require("knitr")
knitr::opts_chunk$set(warning=FALSE, message=FALSE, fig.align='center')

library(ggplot2)
library(broom.mixed)
library(broom)
library(tidyverse)
library(tidymodels)
library(plyr)
library(dplyr)
library(plotrix)
library(lubridate)
library(effects)
library(RCurl)
library(chron)
library(lubridate)
library(lme4)
library(emmeans)
library(multcomp)
library(devtools)
library(logihist)
library(popbio)
library(car)
library(lmerTest)
library(cowplot)
library(sjPlot)
library(sjmisc)
library(RgoogleMaps)
library(SDMTools)
library(rgdal)
library(maps)
library(png)
library(ecodist)
library(gridExtra)

Project background

Montipora capitata was sampled at four reef locations in Kāne‘ohe Bay, O‘ahu, spanning a latitudinal (north-south) gradient of oceanic input and seawater residence and nearshore proximity (east-west)

We measured changes in total biomass (ash-free dry weight), symbiont densiteis, photopigments (total chlorophylls (a+c2)), and stable isotopes (δ13C, δ15N) of the coral host, its endosymbionts Symbiodinium, and the coral skeleton.

Site map

# attach map GPS data
HI<-readOGR("data/coast_n83.shp", "coast_n83") # in meters
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83"))

sites<-read.csv("data/environmental/Reefs_lat_long.csv", header=TRUE, na.string=NA)
# Make figure and export

###############

par(mar=c(1,1,1,1), pty="sq")
plot(HI, xlim=c(-157.86, -157.75), ylim=c(21.42, 21.5), lwd=1, col="gray87")

# adding tick marks and axis lat-long
axis(side=1, at=seq(-157.9, -157.70, 0.05), line=0, tck=0.02, labels=FALSE) # bottom
text(y=21.409, x=seq(-157.9, -157.70, 0.05), labels=seq(-157.9, -157.70, 0.05), cex=0.4)

axis(side=3, at=seq(-157.9, -157.70, 0.05), line=0, tck=0.02, labels=FALSE) # top

axis(side=2, at=seq(21.35, 21.55, 0.05), line=0, tck=0.02, labels=FALSE) # right
text(x=-157.857, y=seq(21.35, 21.55, 0.05), labels=seq(21.35, 21.55, 0.05),cex=0.4)

axis(side=4, at=seq(21.35, 21.55, 0.05), line=0, tck=0.02, labels=FALSE) # left

## add points to map, needs to be 'longitude-latitude formatted'
points(sites[,c(3,2)], pch=21, cex=1.3, col="black", bg="salmon")
text(sites[,c(3,2)], labels=c("SE\npatch reef", "SW\nfringe reef", "NE\npatch reef", "NW\nfringe reef"), pos=c(3,3,1,3), cex=0.5)

par(new=T, mar=c(15,15,1,1))
plot(HI, xlim=c(-158.27, -157.6), ylim=c(21.25, 21.72), lwd=0.4, col="gray", bg="white")  # Oahu
rect(-157.87, 21.39, -157.71, 21.53, lwd=0.8)
box()
**Figure 1a.** *Study site locations in Kāne‘ohe Bay on the windward side of the island of O'ahu, Hawai'i*

Figure 1a. Study site locations in Kāne‘ohe Bay on the windward side of the island of O’ahu, Hawai’i 

##### save the figure and export to directory? ####
dev.copy(pdf, "figures/environmental/KBaymap.pdf", height=4.6, width=5)
dev.off

Figure Montipora capitata in Kāne’ohe Bay, O’ahu, Hawai’i. PC: C. Wall

Environmental

Sea level correction

Corals were collected across in summer and winter 2016 across several different days with each period. To correct for differences in the depth of coral colonies we used NOAA tide data to correct all depth measurements to Mean Sea Level (MSL). This approach was done in Innis et al. 2018 on a larger collection of the coral samples–which have been subset in the current study.

#########################
# Sea level correction
#########################
rm(list=ls())

#### attach data files
data<-read.csv("data/mastersheet_PanKBAY.csv") # master file
qPCR.Innis<-read.csv("data/PanKBay_summer_qPCR.csv") # qPCR from summer (Innis et al. 2018)

JuneTide=read.csv("data/environmental/sea level/Station_1612480_tide_ht_20160601-20160630.csv")
JulyTide=read.csv("data/environmental/sea level/Station_1612480_tide_ht_20160701-20160731.csv")
AugustTide=read.csv("data/environmental/sea level/Station_1612480_tide_ht_20160801-20160831.csv")
DecemberTide=read.csv("data/environmental/sea level/Station_1612480_tide_ht_20161201-20161222.csv")

Tide<-rbind(JuneTide, JulyTide, AugustTide, DecemberTide) # all MSL in meters

Tide$Time <- as.POSIXct(Tide$TimeUTC, format="%Y-%m-%d %H:%M:%S", tz="UTC")
attributes(Tide$Time)$tzone <- "Pacific/Honolulu"
# use attr(as.POSIXlt(Tide$Time),"tzone") to confirm in PST

Tide<-Tide[, c(-5:-8)] #remove unnecessary columns

data$date.time<- as.POSIXct(paste(as.character(data$Date), as.character(data$Time.of.collection)),
                 format="%m/%d/%y %H:%M", tz="Pacific/Honolulu") # make sure time in HST for master

data$Time.r<-round_date(data$date.time,unit="6 minutes")
Tide$Time.r <- Tide$Time
data<-merge(data, Tide, by="Time.r", all.x=T)
data$newDepth <- data$Depth..m-data$TideHT # new depth in METERS

sum.data<- data[(data$Season=="summer"),]; sum.depth.dif<- sum.data$newDepth-sum.data$Depth..m
win.data<- data[(data$Season=="winter"),]; win.depth.dif<- win.data$newDepth-win.data$Depth..m

#### make a plot of differences
# * note that plot works in Knit and alone in code chunk but isn't running if you "run all"*

plot(sum.depth.dif~sum.data$Depth..m, col="pink4", pch=21, bg="plum", cex.axis=0.8, cex.lab=0.8, ylim=c(-0.3, 0.3), xlim=c(0,10),
     xlab="Field depth observation (m)", 
     ylab="Corrected depth difference (m)")
abline(h=0, lty=2, col="gray")

par(new=T)
plot(win.depth.dif~win.data$Depth..m, col="deepskyblue4", pch=21, bg="deepskyblue1", 
     ylim=c(-0.3, 0.3), xlim=c(0,10), xaxt="n", yaxt="n", xlab="", ylab="")
abline(h=0, lty=2, col="gray")
legend("topright", c("Summer", "Winter"), pch=c(21,21), pt.bg=c("plum", "deepskyblue1"), col=c("pink4", "deepskyblue4"), cex=0.9, y.intersp = 1, x.intersp = 0.6, bty="n", inset=c(0, 0))
**Figure** *The differences in the depth where corals were collected in the field (x-axis) relative to the difference in the corrected depth (y-axis) using Mean Sea Level at the time of collection. Tidal correction from NOAA tide data.*

Figure The differences in the depth where corals were collected in the field (x-axis) relative to the difference in the corrected depth (y-axis) using Mean Sea Level at the time of collection. Tidal correction from NOAA tide data.

Light parameters

DLI 2m all sites

This figure shows the daily integrated light (DLI mol photon m-2 d-1) at long-term light loggers at 2m depth near the locations where corals were collected in northern and southern sectors of Kāne’ohe Bay near shore (east) further off shore (west).

# DLI
##### all DLI for graphs
files <- list.files(path="data/environmental/temp and light/Jun_DecPAR/all DLI", pattern = "csv$", full.names = T)
tables <- lapply(files, read.csv, header = TRUE)
DLI<-do.call(rbind, tables)
DLI=DLI[-1] # remove junk column
DLI$Date<-as.Date(DLI$Date) # fix date
DLI<-DLI[order(DLI$Date),] # order by Date

# make a date sequence for entire study
all.date<-as.data.frame(seq(as.Date("2016-06-10"), as.Date("2017-01-12"), "days"))
colnames(all.date)="Date"

# merge the DLI data and the date sequence to make a complete df through time
DLI.3site<-merge(all.date, DLI, by="Date", all.x=T)
R10<-read.csv("data/environmental/temp and light/Jun_DecPAR/all DLI/Reef 10/DLI.Rf102016.csv")
R10<-R10[-1]; R10$Date<-as.Date(R10$Date)

DLI.4site<-merge(DLI.3site, R10, by="Date", all.x=T)
DLI.4site$month <- months(as.Date(DLI.4site$Date)) # makes a month column
DLI.4site<-DLI.4site[, c(1,6, 2:5)]

# write.csv(DLI.4site, "data/environmental/temp and light/Jun_DecPAR/All.DLI.csv")

##### Figure
All.DLI<-DLI.4site
reefcols=c("mediumseagreen", "dodgerblue", "salmon", "orchid")
par(mar=c(2,3.6,1,1.3), mgp=c(2,0.5,0))
layout(matrix(c(1,2,3,4), 2,2, byrow=TRUE))

plot(Rf44.mid.DLI~Date, type="l", All.DLI, ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), ylim=c(0, 30), xaxt="n", xlab="Time", col=reefcols[1])
legend("topright", lty=1, col=reefcols[1], legend="Northwest", lwd=1.5, bty="n", cex=0.8, 
       x.intersp = 0.3, inset=c(0.02, -0.02))
axis.Date(1, at=seq(min(All.DLI$Date), max(All.DLI$Date), by="2 mo"), format="%b '%y")

plot(Rf42.mid.DLI~Date, type="l", All.DLI, ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), ylim=c(0, 30), xaxt="n", xlab="Time", col=reefcols[2])
legend("topright", lty=1, col=reefcols[2], legend="Northeast", lwd=1.5, bty="n", cex=0.8, 
       x.intersp = 0.3, inset=c(0.02, -0.02))
axis.Date(1, at=seq(min(All.DLI$Date), max(All.DLI$Date), by="2 mo"), format="%b '%y")

plot(Rf10.mid.DLI~Date, type="l", All.DLI, ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1, ")", sep=""))), ylim=c(0, 30), xaxt="n", xlab="Time", col=reefcols[3])
legend("topright", lty=1, col=reefcols[3], legend="Southwest", lwd=1.5, bty="n", cex=0.8, 
       x.intersp = 0.3, inset=c(0.02, -0.02))
axis.Date(1, at=seq(min(All.DLI$Date), max(All.DLI$Date), by="2 mo"), format="%b '%y")

plot(HIMB.mid.DLI~Date, type="l", All.DLI, ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1, ")", sep=""))), ylim=c(0, 30), xaxt="n", xlab="Time", col=reefcols[4])
legend("topright", lty=1, col=reefcols[4], legend="Southeast", lwd=1.5, bty="n", cex=0.8, 
       x.intersp = 0.3, inset=c(0.02, -0.02))
axis.Date(1, at=seq(min(All.DLI$Date), max(All.DLI$Date), by="2 mo"), format="%b '%y")
**Figure** *DLI at the four reef areas from summer to winter 2016 at 2m depths. Data at 2m was recorded from loggers*

Figure DLI at the four reef areas from summer to winter 2016 at 2m depths. Data at 2m was recorded from loggers 

dev.copy(pdf, "figures/environmental/All.DLI.pdf", width=6.5, height=4)
dev.off()

Model results testing differences in DLI at 2 m depth at all sites.

mod.light<-DLI.4site
Rf44DLI.mod<-mod.light[1:3]; Rf44DLI.mod$Location<-"Northwest"; 
    colnames(Rf44DLI.mod)<-c("Date", "month", "DLI", "Location")
Rf42DLI.mod<-mod.light[,c(1:2,4)]; Rf42DLI.mod$Location<-"Northeast"; 
    colnames(Rf42DLI.mod)<-c("Date", "month", "DLI", "Location")
HIMBDLI.mod<-mod.light[,c(1:2,5)]; HIMBDLI.mod$Location<-"Southeast"; 
    colnames(HIMBDLI.mod)<-c("Date", "month", "DLI", "Location")
Rf10DLI.mod<-mod.light[,c(1:2,6)]; Rf10DLI.mod$Location<-"Southwest"; 
    colnames(Rf10DLI.mod)<-c("Date", "month", "DLI", "Location")

All.DLI.mod<-rbind(Rf44DLI.mod, Rf42DLI.mod, HIMBDLI.mod, Rf10DLI.mod)
All.DLI.mod$Location<-as.factor(All.DLI.mod$Location)
All.DLI.mod$month<-as.factor(All.DLI.mod$month)
All.DLI.mod$Season<-ifelse(All.DLI.mod$month== "June" | 
                             All.DLI.mod$month== "July" |
                                All.DLI.mod$month== "August" |
                                  All.DLI.mod$month== "September", "dry season", "wet season")
All.DLI.mod$Season<-as.factor(All.DLI.mod$Season)

mod<-lmer(DLI~Location*Season + (1|Date), data=All.DLI.mod); anova(mod, type=2)
## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Location        4378.8 1459.58     3 529.66 134.673 < 2.2e-16 ***
## Season          1040.9 1040.91     1 210.48  96.043 < 2.2e-16 ***
## Location:Season  490.9  163.62     3 531.35  15.097 1.923e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(mod))
posthoc<-emmeans(mod, ~Location|Season)
CLD(posthoc, Letters=letters)
## Season = dry season:
##  Location     emmean        SE     df  lower.CL  upper.CL .group
##  Southwest  8.149681 0.4014748 504.55  7.360912  8.938449  a    
##  Northwest 10.829854 0.4450350 590.91  9.955811 11.703896   b   
##  Southeast 14.748648 0.4016675 503.81 13.959498 15.537797    c  
##  Northeast 14.997825 0.4601858 615.62 14.094100 15.901549    c  
## 
## Season = wet season:
##  Location     emmean        SE     df  lower.CL  upper.CL .group
##  Southwest  3.901772 0.4473755 539.33  3.022959  4.780584  a    
##  Southeast  8.607241 0.4662620 551.22  7.691373  9.523109   b   
##  Northwest  9.128302 0.4662620 551.22  8.212435 10.044170   b   
##  Northeast  9.448888 0.4707121 559.55  8.524310 10.373467   b   
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
posthoc<-emmeans(mod, ~Season)
CLD(posthoc, Letters=letters)
##  Season        emmean        SE     df  lower.CL  upper.CL .group
##  wet season  7.771551 0.3348801 206.76  7.111333  8.431768  a    
##  dry season 12.181502 0.3071824 206.48 11.575885 12.787118   b   
## 
## Results are averaged over the levels of: Location 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
**Figure** DLI model output

Figure DLI model output

PAR and DLI at 3 depths

With Depth as a factor (<1m, 2m, 8m), we can produce a continuous graph of light (PAR or DLI) at 2m and use model intercepts to predict the light at the shallow and deep zones (ca. <1m and 8m, respectfully). Here, the model is \(lm(log.PAR\) ~ \(Depth.factor)\).

# Plot it!!
# Reef 44
par(mfrow=c(2,2), mar=c(4,5,2,2))

plot(y=All.PAR.df$Rf44.shall, x=All.PAR.df$timestamp, type="l", col="palegreen3", 
     cex.main=0.7, ylab=(expression(paste("PAR"~(mu*mol~photons~m^-2~s^-1), sep="")))
     , xlab="", main="Northwest", cex.main=1, ylim=c(0, 2400))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf44.mid, x=All.PAR.df$timestamp, type="l", col="gray"))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf44.deep, x=All.PAR.df$timestamp, type="l", col="palegreen4"))
legend("topright", lty=1, col=c("palegreen3", "gray", "palegreen4"), legend=c("<1m", "2m", "8m"), 
       lwd=2, bty="n", inset=c(0, -0.01), seg.len=1, cex=0.9, x.intersp=0.5, y.intersp=0.8)


# Plot it!!
# Reef 42
plot(y=All.PAR.df$Rf42.shall, x=All.PAR.df$timestamp, type="l", col="mediumturquoise", 
     cex.main=0.7, ylab=(expression(paste("PAR"~(mu*mol~photons~m^-2~s^-1), sep="")))
     , xlab="", main="Northeast", cex.main=1, ylim=c(0, 2400))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf42.mid, x=All.PAR.df$timestamp, type="l", col="gray"))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf42.deep, x=All.PAR.df$timestamp, type="l", col="dodgerblue2"))
legend("topright", lty=1, col=c("mediumturquoise", "gray", "dodgerblue2"), legend=c("<1m", "2m", "8m"), 
       lwd=2, bty="n", inset=c(0, -0.01), seg.len=1, cex=0.9, x.intersp=0.5, y.intersp=0.8)

# Plot it!!
# Reef 10
plot(y=All.PAR.df$Rf10.shall, x=All.PAR.df$timestamp, type="l", col="orangered", 
     cex.main=0.7, ylab=(expression(paste("PAR"~(mu*mol~photons~m^-2~s^-1), sep="")))
     , xlab="Dates", main="Southwest", cex.main=1, ylim=c(0, 2400))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf10.mid, x=All.PAR.df$timestamp, type="l", col="gray"))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf10.deep, x=All.PAR.df$timestamp, type="l", col="lightsalmon3"))
legend("topright", lty=1, col=c("orangered", "gray", "lightsalmon3"), legend=c("<1m", "2m", "8m"), 
       lwd=2, bty="n", inset=c(0, -0.01), seg.len=1, cex=0.9, x.intersp=0.5, y.intersp=0.8)

# Plot it!
# HIMB
plot(y=All.PAR.df$HIMB.shall, x=All.PAR.df$timestamp, type="l", col="mediumorchid2", 
     cex.main=0.7, ylab=(expression(paste("PAR"~(mu*mol~photons~m^-2~s^-1), sep="")))
     , xlab="Dates", main="Southeast", cex.main=1, ylim=c(0, 2400))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$HIMB.mid, x=All.PAR.df$timestamp, type="l", col="gray"))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$HIMB.deep, x=All.PAR.df$timestamp, type="l", col="plum4"))
legend("topright", lty=1, col=c("mediumorchid2", "gray", "plum4"), legend=c("<1m", "2m", "8m"), 
       lwd=2, bty="n", inset=c(0, -0.01), seg.len=1, cex=0.9, x.intersp=0.5, y.intersp=0.8)
**Figure** *PAR at the four reef areas from summer to winter 2016 at three depths (<1m, 2m, 8m). Data at 2m was recorded from loggers; data at <1m and 8m were modeled using light attenuation calculations*

Figure PAR at the four reef areas from summer to winter 2016 at three depths (<1m, 2m, 8m). Data at 2m was recorded from loggers; data at <1m and 8m were modeled using light attenuation calculations 

dev.copy(pdf, "figures/environmental/PAR.all depths.pdf", height=7, width=8)
dev.off()
#######
####### calculate DLI and compare
#######
df<-All.PAR.df
df$timestamp<-strptime(df$timestamp, format="%Y-%m-%d %H:%M:%S")
df$timestamp<-as.Date(df$timestamp, format="%Y-%m-%d") # change to DATE ONLY format to calulate range, means
df[is.na(df)] <- 0


df.split <- split(df, f=df$timestamp
                  < as.Date("2016-06-10", format="%Y-%m-%d")) # split df by date

df.dli<-aggregate(data.frame(Rf42.shall.DLI=df.split[[1]]$Rf42.shall*0.0864,
                             Rf42.mid.DLI=df.split[[1]]$Rf42.mid*0.0864,
                             Rf42.deep.DLI=df.split[[1]]$Rf42.deep*0.0864,
                             Rf44.shall.DLI=df.split[[1]]$Rf44.shall*0.0864,
                             Rf44.mid.DLI=df.split[[1]]$Rf44.mid*0.0864,
                             Rf44.deep.DLI=df.split[[1]]$Rf44.deep*0.0864,
                             HIMB.shall.DLI=df.split[[1]]$HIMB.shall*0.0864,
                             HIMB.mid.DLI=df.split[[1]]$HIMB.mid*0.0864,
                             HIMB.deep.DLI=df.split[[1]]$HIMB.deep*0.0864,
                             Rf10.shall.DLI=df.split[[1]]$Rf10.shall*0.0864,
                             Rf10.mid.DLI=df.split[[1]]$Rf10.mid*0.0864,
                             Rf10.deep.DLI=df.split[[1]]$Rf10.deep*0.0864),
                  by=list(Date=df.split[[1]]$timestamp), FUN=mean)


# make a date sequence for entire study
all.date.time<-as.data.frame(seq(
  from=as.POSIXct("2016-06-10", tz="HST"),
  to=as.POSIXct("2017-01-12", tz="HST"),
  by="1 d"))  
colnames(all.date.time)[1]<-"Date"
all.date.time$Date<-strptime(all.date.time$Date, format="%Y-%m-%d")
all.date.time$Date<-as.Date(all.date.time$Date, format="%Y-%m-%d")

df.dli<-merge(all.date.time, df.dli, by="Date", all.x=T)
df.dli[df.dli<=0] <- NA


#####################
# Plot it!!
# Reef 44
par(mfrow=c(2,2), mar=c(4,5,2,2))

plot(y=df.dli$Rf44.shall, x=df.dli$Date, type="h", col="palegreen3", 
     cex.main=0.7, ylab=(expression(paste("DLI"~(mol~photons~m^-2~d^-1), sep="")))
     , xlab="", main="Northwest", cex.main=1, ylim=c(0, 50))
par(new=T)
with(df.dli, lines(y=df.dli$Rf44.mid, x=df.dli$Date, type="h", col="gray"))
par(new=T)
with(df.dli, lines(y=df.dli$Rf44.deep, x=df.dli$Date, type="h", col="palegreen4"))
legend("topright", lty=1, col=c("palegreen3", "gray", "palegreen4"), legend=c("<1m", "2m", "8m"), 
       lwd=2, bty="n", inset=c(0.01, 0), seg.len=1, cex=0.9, x.intersp=0.2, y.intersp=1)


# Plot it!!
# Reef 42
plot(y=df.dli$Rf42.shall, x=df.dli$Date, type="h", col="mediumturquoise", 
     cex.main=0.7, ylab=(expression(paste("DLI"~(mol~photons~m^-2~d^-1), sep="")))
     , xlab="", main="Northeast", cex.main=1, ylim=c(0, 50))
par(new=T)
with(df.dli, lines(y=df.dli$Rf42.mid, x=df.dli$Date, type="h", col="gray"))
par(new=T)
with(df.dli, lines(y=df.dli$Rf42.deep, x=df.dli$Date, type="h", col="dodgerblue2"))
legend("topright", lty=1, col=c("mediumturquoise", "gray", "dodgerblue2"), legend=c("<1m", "2m", "8m"), 
       lwd=2, bty="n", inset=c(0.01, 0), seg.len=1, cex=0.9, x.intersp=0.2, y.intersp=1)

# Plot it!!
# Reef 10
plot(y=df.dli$Rf10.shall, x=df.dli$Date, type="h", col="orangered", 
     cex.main=0.7, ylab=(expression(paste("DLI"~(mol~photons~m^-2~d^-1), sep="")))
     , xlab="Dates", main="Southwest", cex.main=1, ylim=c(0, 50))
par(new=T)
with(df.dli, lines(y=df.dli$Rf10.mid, x=df.dli$Date, type="h", col="gray"))
par(new=T)
with(df.dli, lines(y=df.dli$Rf10.deep, x=df.dli$Date, type="h", col="lightsalmon3"))
legend("topright", lty=1, col=c("orangered", "gray", "lightsalmon3"), legend=c("<1m", "2m", "8m"), 
       lwd=2, bty="n", inset=c(0.01, 0), seg.len=1, cex=0.9, x.intersp=0.2, y.intersp=1)

# Plot it!
# HIMB
plot(y=df.dli$HIMB.shall, x=df.dli$Date, type="h", col="mediumorchid2", 
     cex.main=0.7, ylab=(expression(paste("DLI"~(mol~photons~m^-2~d^-1), sep="")))
     , xlab="Dates", main="Southeast", cex.main=1, ylim=c(0, 50))
par(new=T)
with(df.dli, lines(y=df.dli$HIMB.mid, x=df.dli$Date, type="h", col="gray"))
par(new=T)
with(df.dli, lines(y=df.dli$HIMB.deep, x=df.dli$Date, type="h", col="plum4"))
legend("topright", lty=1, col=c("mediumorchid2", "gray", "plum4"), legend=c("<1m", "2m", "8m"), 
       lwd=2, bty="n", inset=c(0.01, 0), seg.len=1, cex=0.9, x.intersp=0.2, y.intersp=1)
**Figure** *DLI at the four reef areas from summer to winter 2016 at three depths (<1m, 2m, 8m). Data at 2m was recorded from loggers; data at <1m and 8m were modeled using light attenuation calculations*

Figure DLI at the four reef areas from summer to winter 2016 at three depths (<1m, 2m, 8m). Data at 2m was recorded from loggers; data at <1m and 8m were modeled using light attenuation calculations 

dev.copy(pdf, "figures/environmental/DLIcalc.all depths.pdf", height=7, width=8)
dev.off()
DLI.summary<-read.csv("data/environmental/temp and light/Jun_DecPAR/DLI.3depth.summary.csv")
DLI.summary$Location<-factor(DLI.summary$Location, levels=c("NW", "NE", "SW", "SE"))
DLI.summary$Depth.zone<-factor(DLI.summary$Depth.zone, levels=c("shallow", "mid", "deep"))

#####
# figure formatting script
Fig.formatting<-(theme_classic()) +
  theme(text=element_text(size=10),
        axis.line=element_blank(),
        legend.justification=c(1,1), legend.position=c(1,1),
        legend.background = element_rect("transparent", colour=NA),
        legend.text=element_text(size=10),
        legend.title = element_blank(),
        panel.border = element_rect(fill=NA, colour = "black", size=0.8),
        aspect.ratio=1, 
        axis.ticks = element_line(colour = 'black', size = 0.4),
        axis.ticks.length=unit(0.25, "cm"),
        axis.text.y=element_text(
          margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=10), 
        axis.text.x=element_text(
          margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=10)) +
theme(legend.spacing.x = unit(0.2, 'cm'),
      legend.key.size = unit(0.4, "cm"))

pd <- position_dodge(0.7) #offset for error bars and columns
#####

Location.label<-c("NW", "NE", "SW", "SE")
Depth.label<-c("<1m", "2m", "8m")
plot_colors2<-c("lightskyblue", "dodgerblue", "deepskyblue4")
plot_stat<-c("slategray", "dodgerblue2", "steelblue1", "slategray", "dodgerblue2", "steelblue1",
             "slategray", "dodgerblue2", "steelblue1","slategray", "dodgerblue2", "steelblue1")

# DLI by site and depth
Fig.DLI<-ggplot(DLI.summary, aes(x=Location, y=DLI.mean, fill=Depth.zone)) +
  geom_bar(stat="identity", col=plot_stat, position = pd, width=0.7, size=0.3) +
  geom_errorbar(aes(ymin=DLI.mean-DLI.se, ymax=DLI.mean+DLI.se),
                size=0.4, width=0, position=pd, col=plot_stat) +
  xlab("Locations") +
  ylab(expression(paste("DLI", ~(mol~photons~m^-2~d^-1), sep=""))) +
  scale_x_discrete(labels=Location.label) +
  scale_y_continuous(expand=c(0,0), limits=c(0, 35)) +
  scale_fill_manual(values=alpha(c(plot_colors2), 0.6), 
                    labels=Depth.label)+ Fig.formatting; Fig.DLI
**Figure** Daily Light Integral (DLI) for four reef locations at three depths. DLI estimated from light attenuation coefficents using depth as a fixed effect against baseline loggers at 2 m depth

Figure Daily Light Integral (DLI) for four reef locations at three depths. DLI estimated from light attenuation coefficents using depth as a fixed effect against baseline loggers at 2 m depth 

ggsave(file="figures/environmental/DLI.bar.pdf", height=4, width=4)

Attenuation and Light-at-Depth

A periodic deployment of loggers (October and November) at ca. <1m and ca. 7m were used to model the relationship between 2m light and light at other depths.

Using the logger data, the difference in light at 1 and 7m was calculated relative to 2m (i.e, PAR baseline (2m) - PAR other), then the difference in depth at 1 and 7m was calculated realtive to logger depth i.e, depth logger other - depth logger (2m baseline). A no-intercept model was used to fit difference in log(DLI) and depth (i.e, log(delta-DLI) and delta-depth) with depth as a continous/numeric variable to determine kd: \(lm(delta.log.DLI\)~\(delta.Depth + 0)\).

Site-specific kd values were used to generate a continuous variable of “light at depth”. At each time period (summer or winter) the daily light integral was averaged over 3 months: summer 2016 (June, July, August) and winter 2016 (November, December, January). Discrete monthly seasonal-mean DLI means was used to generate an estimate of light at depth for corals at each site, as this integrates the site and seasoanl variance.

Light at Depth was calculated following Beer-Lambert: \(Ez{_d} = Ez_{_0}e^{(-k_d *d)}\), where \(d\) is depth, \(k_d\) is the attenuation coefficient for each site, \(Ez{_d}\) is light at depth \(d\) and \(Ez{_0}\) is measured light. Calculations were relative to light at 2m depth (i.e., \(Ez_{_0}\) was the light at 2m depth) and \(d\) was the difference between coral depth relative to depth of the logger (i.e., 2m).

########################## # LIGHT AT DEPTH
#############
###### 
# using new depth and the site-specific kd (light attenuation) determine approximate "light at depth" for each colony sample. 

######

### attach necessary files
logger.depths<-read.csv("data/environmental/temp and light/light.logger.depths.csv") # depths for loggers
kd.all<-read.csv("data/environmental/kd.all.csv") # kds for each site

# Seasonal DLI used for "period of collection" light levels
month.DLI<-read.csv("data/environmental/temp and light/Jun_DecPAR/All.DLI_long.csv")

# corals collected in June-July-August use summer time DLI for these months as indicator of average DLI
summer.DLI<-month.DLI[(month.DLI$Month=="June" | month.DLI$Month=="July" | month.DLI$Month=="August"),]
winter.DLI<-month.DLI[(month.DLI$Month=="November" | month.DLI$Month=="December" |month.DLI$Month=="January"),]

# summer mean and SE dataframe
sum.mean<-aggregate(DLI~Site, summer.DLI, mean)
sum.SE<-aggregate(DLI~Site, summer.DLI, std.error)
sum.light.df<-cbind(sum.mean, sum.SE[2])
colnames(sum.light.df)=c("Site", "mean.DLI", "se.DLI")
sum.light.df$Season<-as.factor("summer")

# winter mean and SE dataframe
wint.mean<-aggregate(DLI~Site, winter.DLI, mean)
wint.SE<-aggregate(DLI~Site, winter.DLI, std.error)
wint.light.df<-cbind(wint.mean, wint.SE[2])
colnames(wint.light.df)=c("Site", "mean.DLI", "se.DLI")
wint.light.df$Season<-as.factor("winter")

season.DLI<-rbind(sum.light.df[,c(4,1,2:3)], wint.light.df[,c(4,1,2:3)]) # compiled means for DLI at ~2m
season.DLI$Location<- as.factor(c("SE", "SW", "NE", "NW", "SE", "SW", "NE", "NW")); season.DLI<-season.DLI[,c(1,5,2:4)]

write.csv(season.DLI, "data/environmental/season.DLI.csv")
#######################################
### make new dataframe for calculations
df.light<- data[, c("Season", "Reef", "Location", "newDepth")]

# make a column of "depth differences" relative to where ~2m logger was deployed
df.light$depth.diff<-ifelse(df.light$Reef=="F1-R46", df.light$newDepth-1.8288,
                          ifelse(df.light$Reef=="F8-R10", df.light$newDepth-1.8288, 
                                 ifelse(df.light$Reef=="HIMB", df.light$newDepth-1.8288, 
                                        df.light$newDepth-2.4384))) # last statement for remaining site (Rf42)

# make a column for sample-specific light at depth (estimate) based on kd
# follow: #with 2m as baseline#  E(depth) = E(2m)*exp(-k_d*(depth - 2m))
# so that:     light at depth = DLI at mid.depth * exp (-kd *(delta shallow-deep in m))


df.light$Light<-ifelse(df.light$Reef=="HIMB" & df.light$Season=="summer", 
                          season.DLI$mean.DLI[1]*exp(-kd.all$kd[1]*df.light$depth.diff), # summer HIMB
                   ifelse(df.light$Reef=="F8-R10" & df.light$Season=="summer", 
                          season.DLI$mean.DLI[2]*exp(-kd.all$kd[2]*df.light$depth.diff), # summer R10
                   ifelse(df.light$Reef=="R42" & df.light$Season=="summer", 
                          season.DLI$mean.DLI[3]*exp(-kd.all$kd[3]*df.light$depth.diff), # summer R42
                   ifelse(df.light$Reef=="F1-R46" & df.light$Season=="summer", 
                                 season.DLI$mean.DLI[4]*exp(-kd.all$kd[4]*df.light$depth.diff), # summer R46
                          
                   ifelse(df.light$Reef=="HIMB" & df.light$Season=="winter", 
                          season.DLI$mean.DLI[5]*exp(-kd.all$kd[1]*df.light$depth.diff), # winter HIMB
                   ifelse(df.light$Reef=="F8-R10" & df.light$Season=="winter", 
                          season.DLI$mean.DLI[6]*exp(-kd.all$kd[2]*df.light$depth.diff), # winter R10
                   ifelse(df.light$Reef=="R42" & df.light$Season=="winter", 
                          season.DLI$mean.DLI[7]*exp(-kd.all$kd[3]*df.light$depth.diff), # winter R42
                          season.DLI$mean.DLI[8]*exp(-kd.all$kd[4]*df.light$depth.diff)) # winter R46
                            ))))))

###### plot of light x depth by season
df.light$Reef <- factor(df.light$Reef, levels = c("F1-R46", "R42", "F8-R10", "HIMB"))
df.light$Location <- factor(df.light$Location, levels=c("NW", "NE", "SW", "SE"))

plot.by.sites=ggplot(df.light, aes(Light)) + geom_density(aes(fill=Location), alpha=0.3, position = 'stack') + scale_x_continuous(limits=c(0, 40)) + ggtitle("Light at Depth by Seasons") 

plot.by.site=ggplot(df.light, aes(Light)) + geom_density(aes(fill=Season), alpha=0.3, position = 'stack') +
  scale_x_continuous(limits=c(0, 40)) + ggtitle("Light at Depth by Seasons") + facet_wrap(~Location, scales="free") + scale_fill_manual(values=c("darkorange1", "dodgerblue1"))

######
# can loop as a list
p=ggplot(df.light, aes(Light, fill=Season)) + geom_density(alpha=0.3, position = 'stack') + scale_x_continuous(limits=c(0, 40)) + ggtitle("Light at Depth by Seasons") 

plots=dlply(df.light, .(Location), function(x) p %+% x + facet_wrap(~Location))


###### plot of light x depth by season with exponential curve fitting
Sum<-df.light[(df.light$Season=="summer"),]
Win<-df.light[(df.light$Season=="winter"),]

plot(Light~newDepth, Sum, col="coral", pch=16, xlab="Depth (m)", ylab="DLI at coral", 
     ylim=c(0, 30), xlim=c(0, 10))
summod<-lm(log(Light)~newDepth, Sum)
curve(exp(coef(summod)["(Intercept)"]+coef(summod)["newDepth"]*x), add=TRUE, col="coral", lwd=2)
par(new=T)
plot(Light~newDepth, Win, col="lightblue2", pch=16, xaxt="n", yaxt="n", 
     xlab="", ylab="", ylim=c(0, 30), xlim=c(0, 10))
wintmod<-lm(log(Light)~newDepth, Win)
curve(exp(coef(wintmod)["(Intercept)"]+coef(wintmod)["newDepth"]*x), add=TRUE, col="lightblue2", lwd=2)
legend("topright", legend=c("summer", "winter"), col=c("coral", "lightblue2"), lty=1, lwd=1, pch=16, bty="n")
**Figure** *Relationship between mean seasonal light availability (mean daily light integral) at depth of coral fragment collection*

Figure Relationship between mean seasonal light availability (mean daily light integral) at depth of coral fragment collection

Dissolved nutrients

Dissolved inorganic nutrients in seawater were quantified on water samples (ca. 100 ml) taken from each site and filtered through a GF/F filter and stored in acid-washed bottles and frozen at -20 °C until analyzes. Dissolved inorganic nutrients—ammonium (NH4+), nitrate + nitrite (NO3- + NO2-), phosphate (PO43-), and silicate (Si(OH)4)—were measured at University of Hawai‘i at Mānoa SOEST Laboratory for Analytical Biogeochemistry using a Seal Analytical AA3 HR nutrient autoanalyzer and expressed as μmol l-1.

Silicate showed no effects (p > 0.323). Phosphate increased in the winter (p=0.046) most notably at northern sites. N+N was higher in northern sites relative to southern sites and increased during the winter at all sites compared to the summer. Ammonium increased in the winter as well (p <0.001).

#########################
#########################
# nutrients
nutr<-read.csv("data/environmental/PanKBay_nutrients.csv")

#models
mod<-lm(silicate~Location+Date, data=nutr); Anova(mod, type=2)
mod<-lm(phosphate~Location+Date, data=nutr); Anova(mod, type=2)
mod<-lm(N.N~Location+Date, data=nutr); Anova(mod, type=2)
mod<-lm(ammonium~Location+Date, data=nutr); Anova(mod, type=2)

nutr$Date<-mdy(nutr$Date) # corrects date format
dates<-cbind("10-Aug '16", "20-Dec '16")

RfHIMB<-nutr[(nutr$Reef=="HIMB"), ]
Rf42<-nutr[(nutr$Reef=="R42"), ]
Rf46<-nutr[(nutr$Reef=="F1-46"), ]
RF10<-nutr[(nutr$Reef=="F8-R10"), ]

Reefs<-c("Reef 46", "Reef 42", "Reef 10", "HIMB")
Locations<-c("NW", "NE", "SW", "SE")
plot_colors<-c("mediumseagreen", "dodgerblue", "salmon", "orchid")

###############
# Phosphate
# par(mfrow=c(1,1), mar=c(2,4,1,1), mgp=c(2,0.5,0))

# figure layout
layout(matrix(c(1,2,3,4), nrow=1, byrow=TRUE))
par(mar=c(5,4.5,7,1.5))

###############
# Silicate
plot(y=Rf46$silicate, x=Rf46$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("silicate"~(mu*mol~L^-1), sep="")), ylim=c(0, 15), pch=19, lty=2, cex=1, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis(side=1, at=Rf46$Date, labels=dates, cex.axis=0.8) # plots Reef 46 / NW
#plots Reef 42/ NE
with(Rf46, lines(Rf42$silicate, x=Rf42$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[2]))
#plots Reef 10/ SW
with(Rf46, lines(RF10$silicate, x=RF10$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[3]))
#plots HIMB/ SE
with(RfHIMB, lines(RfHIMB$silicate, x=RfHIMB$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[4]))

###############
# Phosphate
plot(y=Rf46$phosphate, x=Rf46$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("phosphate"~(mu*mol~L^-1), sep="")), ylim=c(0, 0.25), pch=19, lty=2, cex=1, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis(side=1, at=Rf46$Date, labels=dates, cex.axis=0.8) # plots Rf46 /NW
#plots Reef 42 / NE
with(Rf46, lines(Rf42$phosphate, x=Rf42$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[2]))
#plots Reef 10 /SW
with(Rf46, lines(RF10$phosphate, x=RF10$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[3]))
#plots HIMB / SE
with(RfHIMB, lines(RfHIMB$phosphate, x=RfHIMB$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[4]))

###############
# N+N
plot(y=Rf46$N.N, x=Rf46$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("nitrate+nitrite"~(mu*mol~L^-1), sep="")), ylim=c(0, 1.5), pch=19, lty=2, cex=1, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis(side=1, at=Rf46$Date, labels=dates, cex.axis=0.8) # plots Reef 46/ NW
#plots Reef 42/ NE
with(Rf46, lines(Rf42$N.N, x=Rf42$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[2]))
#plots Reef 10 / SE
with(Rf46, lines(RF10$N.N, x=RF10$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[3]))
#plots HIMB / SW
with(RfHIMB, lines(RfHIMB$N.N, x=RfHIMB$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[4]))


###############
# Ammonium
plot(y=Rf46$ammonium, x=Rf46$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("ammonium"~(mu*mol~L^-1), sep="")), ylim=c(0, 1.5), pch=19, lty=2, cex=1, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis(side=1, at=Rf46$Date, labels=dates, cex.axis=0.8) # plots Reef 46/ NW
#plots Reef 42 /NE
with(Rf46, lines(Rf42$ammonium, x=Rf42$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[2]))
#plots Reef 10/ SW
with(Rf46, lines(RF10$ammonium, x=RF10$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[3]))
#plots HIMB /SE
with(RfHIMB, lines(RfHIMB$ammonium, x=RfHIMB$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[4]))
legend("topleft", legend=Locations, col=plot_colors[1:4], lwd=1, pch=19, lty=2, cex=0.8, bty="n", x.intersp=0.2, y.intersp=1.2, inset=c(0.01, 0), seg.len=1.5)

dev.copy(pdf, "figures/environmental/all.nutrients.pdf", encod="MacRoman", height=3, width=7)
dev.off()

Suspended particulates

Gravimetric SPM

Total (organic + inorganic) suspended particulate matter (SPM) was quantified by collecting and filtering 0.5 L of seawater from each location during each season. Seawater was sieved at 243 μm to remove large debris and filtered onto a pre-combusted (4h, 450°C) GF/F filter (0.7 μm). Salts were removed by light rinsing with ddH2O. Filters were dried at 60°C for 24h and re-weighed to closest 0.001 g to obtain total SPM. Masses were normalized to a liter of water.

###### SPM
spm<-read.csv("data/environmental/PanKBay_SPM.csv")
spm$SPM..g.l<-spm$SPM..g.l*1000 # change to mg
spm$Date<-factor(spm$Date, levels=c("8/10/16", "12/20/16"))
colnames(spm)[4]="SPM.mg"
spm$Location<-ordered(spm$Location, c("NW", "NE", "SW", "SE"))

mod<-lm(SPM.mg~Location+Date, data=spm); Anova(mod, type=2)
## Anova Table (Type II tests)
## 
## Response: SPM.mg
##           Sum Sq Df F value    Pr(>F)    
## Location  8.8382  3 10.3062 0.0003025 ***
## Date      1.4114  1  4.9373 0.0386252 *  
## Residuals 5.4312 19                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-emmeans(mod, ~Location)
CLD(posthoc, Letters=letters)
##  Location    emmean       SE df  lower.CL upper.CL .group
##  NE       0.7966668 0.218271 19 0.3398203 1.253513  a    
##  NW       1.6450002 0.218271 19 1.1881536 2.101847  ab   
##  SE       2.2666668 0.218271 19 1.8098203 2.723513   b   
##  SW       2.2850002 0.218271 19 1.8281536 2.741847   b   
## 
## Results are averaged over the levels of: Date 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
plot(allEffects(mod), par.strip.text=list(cex=0.7))

rm(mean) # this will remove 'mean' from library and let the function below proceed.
mean<-aggregate(SPM.mg~Date + Location+Reef, spm, mean)
se<-aggregate(SPM.mg~Date+Location+Reef, spm, std.error)

spm.df<-cbind(mean, se[c(0,4)])
colnames(spm.df)=c("Date", "Location", "Reef", "SPM", "se")
spm.df$Date<-ordered(spm.df$Date, c("8/10/16", "12/20/16"))
spm.df$Reef<-ordered(spm.df$Reef, c("F1-46", "R42", "F8-10", "HIMB"))
spm.df$Location<-ordered(spm.df$Location, c("NW", "NE", "SW", "SE"))


Reefs<-c("Reef 46", "Reef 42", "Reef 10", "HIMB")
Locations<-c("NW", "NE", "SW", "SE")
Date<-c("10-Aug '16", "20-Dec '16")
plot_colors2<-c("mediumturquoise", "skyblue3", "lightcoral", "plum")

#####
# figure formatting script
Fig.formatting<-(theme_classic()) +
  theme(text=element_text(size=10),
        axis.line=element_blank(),
        legend.justification=c(1,1), legend.position=c(1,1),
        legend.background = element_rect("transparent", colour=NA),
        legend.text=element_text(size=10),
        legend.title = element_blank(),
        panel.border = element_rect(fill=NA, colour = "black", size=0.8),
        aspect.ratio=1, 
        axis.ticks = element_line(colour = 'black', size = 0.4),
        axis.ticks.length=unit(0.25, "cm"),
        axis.text.y=element_text(
          margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=10), 
        axis.text.x=element_text(
          margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=10)) +
theme(legend.spacing.x = unit(0.2, 'cm'),
      legend.key.size = unit(0.4, "cm"))

pd <- position_dodge(0.7) #offset for error bars and columns
#####

# SPM (mg/l)
Fig.SPM<-ggplot(spm.df, aes(x=Date, y=SPM, fill=Location)) +
  geom_bar(colour=c("plum","lightcoral", "skyblue3","mediumturquoise","plum","lightcoral", "skyblue3", "mediumturquoise"),
           stat="identity", position = pd, width=0.7, size=0.3) +
  geom_errorbar(aes(ymin=SPM-se, ymax=SPM+se), 
                colour=c("plum","lightcoral", "skyblue3","mediumturquoise","plum","lightcoral", "skyblue3","mediumturquoise"),
                size=0.4, width=0, position=pd) +
  xlab("Sampling points") +
  ylab(expression(paste("SPM", ~(mg~L^-1), sep=""))) +
  scale_x_discrete(labels=Date) +
  scale_y_continuous(expand=c(0,0), limits=c(0, 5)) +
  scale_fill_manual(values=alpha(c(plot_colors2), 0.5), 
                    labels=Locations)+ Fig.formatting; Fig.SPM

ggsave(file="figures/environmental/SPM.pdf", height=4, width=4)

Isotopes of SPM/plankton

Plankton tows and seawater collections to parameterized particulates and plankters as potential heterotrophic end members available to reef corals. Sampling was done at 4 sites where corals were collected [North: (Reef 42, fringe-Reef 46), and South: (HIMB, fringe-Reef 10)], as well as two reefs in central region of the bay where corals were not collected (Central: Reef 25, fringe-Reef 25) to increase spatial resolution of suspended particulate isotope sample sizes. Using size-fractioned materials collected in seawater and plankton tows, we examined the spatial (among reef sites) and temporal patterns (among seasons) in stable isotope values of size-fractioned end members. Carbon and nitrogen isotope values among the 6 reef sites were not significantly different (p > 0.140), season had marginal effects (p > 0.049), whereas fraction influenced isotope values significantly (p< 0.001). Therefore, isotope values were pooled among reefs and seasons to best interpret size-fraction isotope values. Samples were not acidified prior to analysis as to not affect nitrogen isotope values.

SWiso<-read.csv("data/isotopes_SW_all times.csv")

mod<-(lm(d13C~Location+Season+SW.fraction..um, data=SWiso)); Anova(mod, type=2) # no location effect
## Anova Table (Type II tests)
## 
## Response: d13C
##                  Sum Sq Df F value    Pr(>F)    
## Location         17.926  5  1.3418 0.2626013    
## Season            7.921  1  2.9645 0.0914188 .  
## SW.fraction..um  76.419  4  7.1504 0.0001286 ***
## Residuals       130.920 49                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod<-(lm(d15N~Location+Season+SW.fraction..um, data=SWiso)); Anova(mod, type=2) # no location effect
## Anova Table (Type II tests)
## 
## Response: d15N
##                  Sum Sq Df F value    Pr(>F)    
## Location         2.3255  5  1.7291   0.14557    
## Season           1.0935  1  4.0653   0.04927 *  
## SW.fraction..um 30.3773  4 28.2335 3.454e-12 ***
## Residuals       13.1802 49                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-emmeans(mod, ~Season)
CLD(posthoc, Letters=letters)
##  Season emmean         SE df lower.CL upper.CL .group
##  summer   6.25 0.09468949 49 6.059714 6.440286  a    
##  winter   6.52 0.09468949 49 6.329714 6.710286   b   
## 
## Results are averaged over the levels of: Location, SW.fraction..um 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
SWiso$Reef<-factor(SWiso$Reef, levels=c("F1-46", "R42", "F2-R25", "R25", "F8-R10", "HIMB"))
SWiso$Location<-factor(SWiso$Location, levels=c("NW", "NE", "CW", "CE", "SW", "SE"))
SWiso$SW.fraction..um<-factor(SWiso$SW.fraction..um, levels=c(">243", "<243", "100-243", "10-100", "0-10"))

winter.data<-SWiso[(SWiso$Season=="winter"),]
summer.data<-SWiso[(SWiso$Season=="summer"),]

These box plots show the seasonal isotope values for carbon and nitrogen according to size fractions.

# plots of SW isotopes by Season-- first showing by size, then by site
op<-par(mfrow = c(2,2), mar=c(5,5,2,1))

######### Sizes
# Summer size d13C
plot(d13C~SW.fraction..um, data=summer.data, ylim=c(-25,-10), 
     main=expression(paste("summer"~ delta^{13}, "C")), ylab=expression(paste(delta^{13}, "C (‰, V-PDB)")),
     col="paleturquoise3", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab=expression(paste("Size Fraction (", mu,m,")")))

# Winter size d13C
plot(d13C~SW.fraction..um, data=winter.data, ylim=c(-25,-10), 
     main=expression(paste("winter"~ delta^{13}, "C")), ylab=expression(paste(delta^{13}, "C (‰, V-PDB)")),
     col="paleturquoise3", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab=expression(paste("Size Fraction (", mu,m,")")))

# Summer size d15N
plot(d15N~SW.fraction..um, data=summer.data, ylim=c(3,10), 
     main=expression(paste("summer"~ delta^{15}, "N")), col="plum", cex.axis=0.7, cex.main=1, cex.lab= 0.8,
     xlab=expression(paste("Size Fraction (", mu,m,")")),
     ylab=expression(paste(delta^{15}, "N (‰, air)")))

# Winter size d15N
plot(d15N~SW.fraction..um, data=winter.data, ylim=c(3,10), 
     main=expression(paste("winter"~ delta^{15}, "N")), ylab=expression(paste(delta^{15}, "N (‰, air)")),
     col="plum", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab=expression(paste("Size Fraction (", mu,m,")")))
**Figure** *Size-fractioned sample isotope values stratified by seasons (summer and winter)*

Figure Size-fractioned sample isotope values stratified by seasons (summer and winter)

These figures show the same size-fractioned isotope data as above, but stratified by location and pooled across seasons.

######## Sites
op<-par(mfrow = c(2,2), mar=c(5,5,2,1))

# Summer site d13C
plot(d13C~Location, data=summer.data, ylim=c(-25,-10), 
     main=expression(paste("summer"~ delta^{13}, "C")), ylab=expression(paste(delta^{13}, "C (‰, V-PDB)")),
     col="lightskyblue", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab="Reef Sites")

# Winter site d13C
plot(d13C~Location, data=winter.data, ylim=c(-25,-10), 
     main=expression(paste("winter"~ delta^{13}, "C")), ylab=expression(paste(delta^{13}, "C (‰, V-PDB)")),
     col="lightskyblue", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab="Reef Sites")

# Summer site d15N
plot(d15N~Location, data=summer.data, ylim=c(3,10), 
     main=expression(paste("summer"~ delta^{15}, "N")), col="salmon", cex.axis=0.7, cex.main=1, cex.lab= 0.8,
     ylab=expression(paste(delta^{15}, "N (‰, air)")), xlab="Reef Sites")

# Winter site d15N
plot(d15N~Location, data=winter.data, ylim=c(3,10), 
     main=expression(paste("winter"~ delta^{15}, "N")), ylab=expression(paste(delta^{15}, "N (‰, air)")),
     col="salmon", cex.axis=0.7, cex.main=1, xlab="Reef Sites")
**Figure** *Isotope values in seawater particles stratified by sites (northwest to southeast)*

Figure Isotope values in seawater particles stratified by sites (northwest to southeast)

Pooling the isotope data across sites and seasons (where effects were absent or negligible), this planktonic foodweb for carbon and nitrogen illustrates differences in the food sources and their isotopic composition.

mod<-lm(d13C~SW.fraction..um, data=SWiso); Anova(mod, type=2)
posthoc<-emmeans(mod, ~SW.fraction..um)
CLD(posthoc, Letters=letters)

mod<-lm(d15N~SW.fraction..um, data=SWiso); Anova(mod, type=2)
posthoc<-emmeans(mod, ~SW.fraction..um)
CLD(posthoc, Letters=letters)


rm(mean)
####
#### making scatter for d15N and d13C, pooled across seasons and sites
mix.N.mean<-aggregate(d15N~SW.fraction..um, data=SWiso, mean)
mix.N.se<-aggregate(d15N~SW.fraction..um, data=SWiso, std.error)
mix.C.mean<-aggregate(d13C~SW.fraction..um, data=SWiso, mean)
mix.C.se<-aggregate(d13C~SW.fraction..um, data=SWiso, std.error)
mix.data<-cbind(mix.N.mean, mix.C.mean[c(2,0)], mix.N.se[c(2,0)], mix.C.se[c(2,0)]); colnames(mix.data)=c("fraction", "d15N", "d13C", "d15N.se", "d13C.se")

colors=c("#FF6A6A", "#00B2EE", "#FFB90F", "#3CB371", "#8B7500")
op<-par(mfrow = c(1,1), mar=c(5,4,1,5),xpd=TRUE, pty="sq")

size.labels=c(expression(paste("> 243"~mu,m)), expression(paste("< 243"~mu,m)), expression(paste("100-243"~mu,m)), expression(paste("10-100"~mu,m)), expression(paste("< 10"~mu,m)))

#### make the plot
plot(d15N~d13C, data=mix.data, type="n", ylim=c(4.5,8), xlim=c(-23,-17), tck=-0.03, cex.axis=0.7, cex.lab=0.7,
     ylab=expression(paste(delta^{15}, "N (‰, air)")), xlab=expression(paste(delta^{13}, "C (‰, V-PDB)")))
legend("topleft", inset=c(0.05,0.0), legend=size.labels, col=colors, pch=19, cex=0.6, bty="n", x.intersp=1.8, y.intersp = 1.4)
arrows(mix.data$d13C-mix.data$d13C.se, mix.data$d15N, mix.data$d13C+mix.data$d13C.se, mix.data$d15N, length=0.0, angle=90, code=3, lwd=0.7, col=colors)
arrows(mix.data$d13C, mix.data$d15N-mix.data$d15N.se, mix.data$d13C, mix.data$d15N+mix.data$d15N.se, length=0.0, angle=90, code=3, lwd=0.7, col=colors)
points(d15N~d13C, data=mix.data, pch=19, cex=0.8, col=colors)
**Figure** *Isotope values for size-fractioned seawater particles and plankters pooled across seasons and reef sites*

Figure Isotope values for size-fractioned seawater particles and plankters pooled across seasons and reef sites 

dev.copy(pdf, "figures/environmental/iso.sources.KBay.pdf", encod="MacRoman", height=4, width=4)
dev.off()

Biology

################################
### Biological responses
################################

#data : this is the master file

# add in light at depth column from df.light dataframe
data$Light<-df.light$Light

##### produce a categorical depth bin ####
depth<-data$newDepth
data$depth.bin<-factor(ifelse(depth<2, "<2m", ifelse(depth >2 & depth <4, "2-4m", ifelse(depth >4 & depth <6, "4-6m", ">6m"))), levels=c("<2m", "2-4m", "4-6m", ">6m"))

aggregate(Sample.ID~depth.bin+Season+Location, data, length)
data$depth.bin.small<-factor(ifelse(depth<4, "<4m", ">4m"), levels= c("<4m", ">4m"))

################################################
# calculate, normalized dependent variables
################################################
str(data)
data$cells.ml<-as.numeric(data$cells.ml)

# helpful shorthand
SA<-data$surface.area # surface area in cm2
blastate<-data$total.blastate.ml # tissue slurry blastate in ml

# AFDW.mg. == convert AFDW g to mg, mutiply by blastate volume, divide by cm2
data$biomass<- (data$mg.biomass.ml*blastate)/SA

# Symbiodinium.cells. == cell.ml * blastate / SA
data$zoox<- (data$cells.ml*blastate)/SA

# total chlorophyll == ug.chl.a.ml * blastate + ug.chl.c2.ml * blastate / SA
data$chltot<-(data$ug.chl.a.ml)+(data$ug.chl.c2.ml)*blastate/SA

# pg.chlorophyll.a..cell + pg.chlorophyll.c2..cell == ug.chltot.ml * 10^6 / cells.ml
data$chlcell<- (data$ug.chl.a.ml*10^6+data$ug.chl.c2.ml*10^6)/data$cells.ml

qPCR

####################
# qPCR
#########

# qPCR

library(devtools)
# Use steponeR to import data and calculate proporation of C and D symbionts
source_url("https://raw.githubusercontent.com/jrcunning/steponeR/master/steponeR.R")
Mcap.plates <- list.files(path="data/qPCR", pattern = "txt$", full.names = T); Mcap.plates
Mcap <- steponeR(files=Mcap.plates, delim="\t",
                 target.ratios=c("C.D"),
                 fluor.norm=list(C=2.26827, D=0),
                 copy.number=list(C=33, D=3),
                 ploidy=list(C=1, D=1), 
                 extract=list(C=0.813, D=0.813))

Mcap <- Mcap$result
head(Mcap)

# remove +/-control
Mcap <- Mcap[grep("+C52", Mcap$Sample.Name, fixed=T, invert = T), ]
Mcap <- Mcap[grep("H2O", Mcap$Sample.Name, fixed=T, invert = T), ]

# to remove any early-amplification CT noise
Mcap$C.CT.mean[which(Mcap$C.CT.mean < 15)] <- 0

#Remove failed samples, i.e., those where either C or D were NOT found in both reps
Mcap$fail <- ifelse(Mcap$C.reps < 2 & Mcap$D.reps < 2, TRUE, FALSE)
fails <- Mcap[Mcap$fail==TRUE, ]
Mcap <- Mcap[which(Mcap$fail==FALSE),]

# replace CT means with 'NA' as zero
Mcap$C.CT.mean[is.na(Mcap$C.CT.mean)] <-0
Mcap$D.CT.mean[is.na(Mcap$D.CT.mean)] <-0

Mcap$C.D[is.na(Mcap$C.D)] <- 1 # sets all infinity (= 100% C) to 1.0

# caluclate proportion C and proprtion D where C and D are both present
Mcap$propC<- Mcap$C.D / (Mcap$C.D + 1)
Mcap$propD<- 1 / (Mcap$C.D + 1)

# where C and D are not cooccuring...
# if C.D = 1 = 100% C, make 'PropC' = 1 and 'PropD' = 0
# if C.D = 0 = 100% D, make 'PropD' = 1 and 'PropC' = 0
Mcap$propC[which(Mcap$C.D==1)] <- 1
Mcap$propD[which(Mcap$propC==1)] <- 0
Mcap$propD[which(Mcap$C.D==0)] <- 1

# calculate FOUR COMMUNITY categories: C, C>D, D>C, D
Mcap$Mix <- factor(ifelse(Mcap$propC > Mcap$propD, ifelse(Mcap$propD!= 0, "CD", "C"), ifelse(Mcap$propD > Mcap$propC, ifelse(Mcap$propC!=0, "DC", "D"), NA)), levels=c("C", "CD", "DC", "D"))

# Identify SINGLE dominant symbiont clade: C or D
Mcap$Dom <- factor(substr(as.character(Mcap$Mix), 1, 1))


######## look for duplicates in dataset by year and type of event (bleach/recover)
Mcap[duplicated(Mcap$Sample.Name), ] ## duplicates

# remove duplicated
Mcap<-Mcap[!(Mcap$Sample.Name=="HIMB_15" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="R10_05" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="R42_06" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="R46_01" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="R46_02" & Mcap$File.Name=="Wall_PanKbay_plate2.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="R46_03" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="HIMB_13" & Mcap$File.Name=="Wall_PanKbay_plate2.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="HIMB_14" & Mcap$File.Name=="Wall_PanKbay_plate2.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="R10_15" & Mcap$File.Name=="Wall_PanKbay_plate3.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="R10_15" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="R42_11" & Mcap$File.Name=="Wall_PanKbay_plate2.txt"),] 

# parse Sample ID and Site from "Site.Name"
Mcap<-cbind(Mcap, colsplit(Mcap$Sample.Name, pattern= "_", c("Reef", "Sample.ID")))
Mcap$Season<-as.factor("winter")
Mcap$Reef<-as.factor(Mcap$Reef)

Mcap$Reef<-revalue(Mcap$Reef, c("R10"="F8-R10", "R46"="F1-R46")) # rename factor levels

# make new factors for bay region and reef type
Mcap$Bay.region <- ifelse(Mcap$Reef=="R42" | Mcap$Reef=="F1-R46", "northern", "southern")
Mcap$Reef.type <- ifelse(Mcap$Reef=="R42" | Mcap$Reef=="HIMB", "patch", "fringe")
Mcap$Location <- ifelse(Mcap$Reef=="F1-R46", "NW", ifelse(Mcap$Reef=="R42", "NE", 
                       ifelse(Mcap$Reef=="F8-R10", "SW", "SE")))

### reorder columns and finish
Mcap<-Mcap[, c(17,15,20,19,18,16,1:14)] # reordered to match masterdata, and finish

### structure winter and summer qPCR dataframes to have same columns and combine dataframe 
qPCR.winter<-Mcap[ , (names(Mcap) %in% 
        c("Season", "Reef", "Location", "Reef.type", "Bay.region", "Sample.ID", "propC", "propD", "Mix", "Dom"))]
qPCR.summer<-qPCR.Innis[ , (names(qPCR.Innis) %in% 
        c("Season", "Reef", "Location", "Reef.type", "Bay.region", "Sample.ID", "propC", "propD", "Mix", "Dom"))] 

# merge qPCR files
qPCR.all<-rbind(qPCR.winter, qPCR.summer)

# add to master data
data.all<-merge(data, qPCR.all, by=c("Season", "Reef", "Location", "Reef.type", "Bay.region", "Sample.ID"), all.x=T)

###### remove columns no longer needed, update "Depth" to be tide-corrected depth )= newDepth
data.trim<-data.all[ , !(names(data.all) %in% c("total.blastate.ml", "Date", "Time.of.collection", "Depth..m", "Time.r", "surface.area.cm2", "cells.ml", "ug.chl.a.ml", "ug.chl.c2.ml", "mg.biomass.ml", "host..mass.mg", "host..ugN", "host..ugC", "symb..mass.mg", "symb..ugN", "symb..ugC", "stationId", "datum", "TimeUTC", "TideHT", "Time"))]

data.trim$symb..C.N[data.trim$symb..C.N>=12.520270]=NA # set this outlier to NA

qPCR figure The distribution of C and D symbiont types is modeled using a logistic regression with dominant symbiont community (C vs. D) as a response and Depth (m), Season, and reef Location as a predictor. A generalized linear model with a binomial distribution and logit link function is used. Model selection (by model AIC values) between a fully crossed model (crossed) and additive model (add) were compared. The additive mode best fit the data; main effects were tested in anova with a Chi-square test.

## qPCR figures

#####
model.data<-data.trim[,c(17, 1:5, 18:25, 7:16, 26:29)]

#########################
#########################
#########################

# Plots with DEPTH as the predictor
# qPCR and symbionts

#########################
###### Plot Dominant Symbiont and Depth (both seasons)
Symb<-model.data
Symb$Reef <- factor(Symb$Reef, levels = c("F1-R46", "R42", "F8-R10", "HIMB")) # set levels
Symb$Location <- factor(Symb$Location, levels = c("NW", "NE", "SW", "SE")) # set levels

Symb$Dominant <- ifelse(Symb$Dom=="D", 1, 0)
add <-glm(Dominant~newDepth+Season+Location, family = "binomial", data = Symb)
crossed <-glm(Dominant~newDepth*Season+Location, family = "binomial", data = Symb)
AIC(add, crossed) # crossed model best
##         df      AIC
## add      6 129.7125
## crossed  7 125.5697
anova(crossed, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Dominant
## 
## Terms added sequentially (first to last)
## 
## 
##                 Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                              119     143.06              
## newDepth         1  21.6786       118     121.38 3.224e-06 ***
## Season           1   0.3020       117     121.08    0.5826    
## Location         3   3.3647       114     117.71    0.3387    
## newDepth:Season  1   6.1428       113     111.57    0.0132 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
results.all<-crossed

#summary(results.all)
plot(allEffects(results.all), par.strip.text=list(cex=0.6))

results.all.mod <-glm(Dominant~newDepth, family = "binomial", data = Symb)
fitted.all <- predict(results.all.mod, newdata = list(newDepth=seq(0,10,0.1)), type = "response")

###### summer only symbionts
sum.dat<-Symb[(Symb$Season=="summer"),]
results.sum=glm(Dominant~newDepth, family = "binomial", data = sum.dat)
anova(results.sum, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Dominant
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                        59      67.48              
## newDepth  1    21.71        58      45.77 3.171e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#summary(results.sum)
fitted.sum <- predict(results.sum, newdata = list(newDepth=seq(0,10,0.1)), type = "response")
# plot(fitted.sum, ylab="proportion D", ylim=c(0,1))

###### winter only symbionts
wint.dat<-Symb[(Symb$Season=="winter"),]
results.win=glm(Dominant~newDepth, family = "binomial", data = wint.dat)
anova(results.win, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Dominant
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                        59     74.920           
## newDepth  1   4.6265        58     70.293  0.03148 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#summary(results.win)
fitted.win <- predict(results.win, newdata = list(newDepth=seq(0,10,0.1)), type = "response")
# plot(fitted.win, ylab="proportion D", ylim=c(0,1))

Logistic regression of the probability of Montipora capitata have a clade D-dominant symbiont community (probability of 1.0) versus one dominated by clade C symbiont (probability of 0). Both seasons show similar patterns, although the probability of corals being clade D-dominant increases in winter months as a function of Depth relative to summer months.

######
######
## Figure of Dominant symbiont clades across seasons
## **Note** where points equal 0.0 is 100% D, where they equal 0 is 100% C
par(mar=c(5,4,3,2))
plot(sum.dat$propD~sum.dat$newDepth, xlab="Depth (m)", ylab = "Proportion of Clade D Symbiont", pch=19, col="salmon", xlim=c(0,10), ylim=c(0.0, 1.0), cex=0.5, cex.lab=0.8)
par(new=T)
plot(wint.dat$propD~wint.dat$newDepth, pch=19, col="lightskyblue", xlim=c(0,10), ylim=c(0.0, 1.0), xaxt="n", xlab="", ylab="", cex=0.5, cex.lab=0.8)
lines(fitted.all ~ seq(0,10,0.1), col="gray30", lwd=1, lty=2)
lines(fitted.sum ~ seq(0,10,0.1), col="coral", lwd=1)
lines(fitted.win ~ seq(0,10,0.1), col="lightskyblue", lwd=1)
legend("topright", pch=c(19,19, NA), lty=c(1,1,2), col=c("coral", "lightskyblue", "gray30"), legend=c("Summer", "Winter", "Combined"), lwd=1, bty="n", x.intersp=1, pt.cex=0.8, y.intersp=2, cex=0.5, inset=c(-0.005, -0.01), seg.len=4)
**Figure** Probability of *Montipora capitata* hosting clade D symbionts across a depth gradient in summer and winter seasons. Lines represent model fit to data by logistic regression for summer alone, winter alone, and combined (summer + winter) data

Figure Probability of Montipora capitata hosting clade D symbionts across a depth gradient in summer and winter seasons. Lines represent model fit to data by logistic regression for summer alone, winter alone, and combined (summer + winter) data 

dev.copy(pdf, "figures/symbionts/Symbionts_by_Season_depth.pdf", encod="MacRoman", height=4, width=4)
dev.off()

Seen another way–Histograms. By separating these data into 2 histograms the probability of clade D- (1.0) and clade C-dominance (0.0) can be visualized. Here it is clear that there is a greater probability of finding corals with clade D dominant communites at shallow depths, whereas in winter there is a slightly greater probability of finding clade D in corals with increasing depth.

#######
#######
par(mfrow=c(1,2))
Dom <- subset(Symb, !is.na(newDepth) & !is.na(Dominant))
Dom.sum<-Dom[(Dom$Season=="summer"),]
Dom.win<-Dom[(Dom$Season=="winter"),]

logi.hist.plot(Dom.sum$newDepth, Dom.sum$Dominant, boxp = FALSE, type = "hist", col="coral", xlabel = "Depth (m)", ylabel = "", ylabel2 = "", main="summer")
mtext(side = 2, text = "Probability of Clade D", line = 3, cex = 1)

logi.hist.plot(Dom.win$newDepth, Dom.win$Dominant, boxp = FALSE, type = "hist", col="lightskyblue", xlabel = "Depth (m)", ylabel = "", ylabel2 = "", main="winter")
mtext(side = 4, text = "Frequency", line = 0.5, cex=1)
mtext(side = 4, text = "C                                             D", line = 0.5, cex = 0.8)
**Figure** Histograms of probability of *Montipora capitata* hosting clade D symbionts across a depth gradient among seasons. Lines represent model fit to data by logistic regression

Figure Histograms of probability of Montipora capitata hosting clade D symbionts across a depth gradient among seasons. Lines represent model fit to data by logistic regression 

dev.copy(pdf, "figures/symbionts/Symb_Season_Light_logistic.pdf", encod="MacRoman", height=5, width=8)
dev.off()
#
#########################
#########################

# this code runs the same as above, except using colony LIGHT-AT-DEPTH as a predictor instead of depth

# Light as a predictor
# qPCR and symbionts

#########################
###### Plot Dominant Symbiont and Depth (both seasons)
Symb<-model.data
Symb$Dominant <- ifelse(Symb$Dom=="D", 1, 0)
add<-glm(Dominant~Light+Season+Location, family = "binomial", data = Symb)
crossed<-glm(Dominant~Light*Season*Location, family = "binomial", data = Symb)
#AIC(add, crossed)

anova(add, test = "Chisq")
#summary(add)

results.all<-add

results.all.mod<-glm(Dominant~Light, family = "binomial", data = Symb)
fitted.all <- predict(results.all.mod, newdata = list(Light=seq(0,25,0.1)), type = "response")

###### summer only symbionts
sum.dat<-Symb[(Symb$Season=="summer"),]
results.sum<-glm(Dominant~Light, family = "binomial", data = sum.dat)
#anova(results.sum, test = "Chisq")
#summary(results.sum)
fitted.sum <- predict(results.sum, newdata = list(Light=seq(0,25,0.1)), type = "response")
# plot(fitted.sum, ylab="proportion D", ylim=c(0,1))

###### winter only symbionts
wint.dat<-Symb[(Symb$Season=="winter"),]
results.win<-glm(Dominant~Light, family = "binomial", data = wint.dat)
#anova(results.win, test = "Chisq")
#summary(results.win)
fitted.win <- predict(results.win, newdata = list(Light=seq(0,25,0.1)), type = "response")
# plot(fitted.win, ylab="proportion D", ylim=c(0,1))

######
######
## Figure of Dominant symbiont clades across seasons
## **Note** where points equal 0.0 is 100% D, where they equal 0 is 100% C
par(mar=c(5,4,3,2))
plot(sum.dat$propD~sum.dat$Light, xlab="Light (DLI)", ylab = "Proportion of Clade D Symbiont", pch=19, col="salmon", xlim=c(0,25), ylim=c(0.0, 1.0), cex=0.5, cex.lab=0.8)
par(new=T)
plot(wint.dat$propD~wint.dat$newDepth, pch=19, col="lightskyblue", xlim=c(0,25), ylim=c(0.0, 1.0), xaxt="n", xlab="", ylab="", cex=0.5)
lines(fitted.all ~ seq(0,25,0.1), col="gray30", lwd=1, lty=2)
lines(fitted.sum ~ seq(0,25,0.1), col="coral", lwd=1)
lines(fitted.win ~ seq(0,25,0.1), col="lightskyblue", lwd=1)
legend("bottomright", pch=c(19,19, NA), lty=c(1,1,2), col=c("coral", "lightskyblue", "gray30"), legend=c("Summer", "Winter", "Combined"), lwd=1, bty="n", x.intersp=1, pt.cex=0.7, y.intersp=2, cex=0.5, inset=c(0.1, 0.15), seg.len=4)

dev.copy(pdf, "figures/symbionts/Symbionts_by_Season_light.pdf", encod="MacRoman", height=4, width=4)
dev.off()


#######
#######
par(mfrow=c(1,2))
Dom <- subset(Symb, !is.na(Light) & !is.na(Dominant))
Dom.sum<-Dom[(Dom$Season=="summer"),]
Dom.win<-Dom[(Dom$Season=="winter"),]

logi.hist.plot(Dom.sum$Light, Dom.sum$Dominant, boxp = FALSE, type = "hist", col="coral", xlabel = "Light-at-depth (DLI)", ylabel = "", ylabel2 = "", main="summer")
mtext(side = 2, text = "Probability of Clade D", line = 3, cex = 1)

logi.hist.plot(Dom.win$Light, Dom.win$Dominant, boxp = FALSE, type = "hist", col="lightskyblue", xlabel = "Light-at-depth (DLI)", ylabel = "", ylabel2 = "", main="winter")
mtext(side = 4, text = "Frequency", line = 0.5, cex=1)
mtext(side = 4, text = "C                                             D", line = 0.5, cex = 0.8)
dev.copy(pdf, "figures/symbionts/Symb_Season_Light_logistic.pdf", encod="MacRoman", height=5, width=8)
dev.off()

Model assumptions

Physiology models

Total Biomass

########## ########## 
######### biomass ---- 
Y<-model.data$biomass
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model
## Data: model.data
## Models:
## add: Y ~ Season + Light + Dom + (1 | Location)
## full: Y ~ Season * Light * Dom + (1 | Location)
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## add   6 846.22 862.95 -417.11   834.22                         
## full 10 848.36 876.23 -414.18   828.36 5.8663      4     0.2094
summary(add)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 828.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.12288 -0.72960 -0.07959  0.55985  3.08964 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept)  5.82    2.412   
##  Residual             60.57    7.783   
## Number of obs: 120, groups:  Location, 4
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   25.62373    2.24200  16.00547  11.429 4.14e-09 ***
## Seasonwinter   0.08287    1.61545 115.55128   0.051    0.959    
## Light         -0.04636    0.17286 108.20061  -0.268    0.789    
## DomD           2.87278    1.75150 115.80076   1.640    0.104    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light 
## Seasonwintr -0.593              
## Light       -0.692  0.471       
## DomD         0.131 -0.258 -0.421
ranova(add)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Y ~ Season + Light + Dom + (1 | Location)
##                npar  logLik    AIC    LRT Df Pr(>Chisq)  
## <none>            6 -414.15 840.29                       
## (1 | Location)    5 -416.26 842.52 4.2221  1     0.0399 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(anova(add, type=2), digits=5)
## Type II Analysis of Variance Table with Satterthwaite's method
##         Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## Season   0.159   0.159     1 115.55  0.0026 0.9592
## Light    4.358   4.358     1 108.20  0.0719 0.7890
## Dom    162.958 162.958     1 115.80  2.6902 0.1037
var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.biomass=var.table
var.biomass$response<-"biomass"


plot(allEffects(add), ylab="biomass", par.strip.text=list(cex=0.7))

Symbiont density

######### symbionts-- 
Y<-model.data$zoox
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+Season:Light + (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model
## Data: model.data
## Models:
## add: Y ~ Season + Light + Dom + Season:Light + (1 | Location)
## full: Y ~ Season * Light * Dom + (1 | Location)
##      Df    AIC    BIC   logLik deviance  Chisq Chi Df Pr(>Chisq)
## add   7 16.414 35.926 -1.20681   2.4136                         
## full 10 21.951 49.826 -0.97577   1.9515 0.4621      3     0.9271
summary(add)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Y ~ Season + Light + Dom + Season:Light + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 30.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.08854 -0.62190 -0.06602  0.57425  2.14357 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.02478  0.1574  
##  Residual             0.05715  0.2391  
## Number of obs: 120, groups:  Location, 4
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)         14.309970   0.101529   6.918170 140.945 3.19e-13 ***
## Seasonwinter        -0.031559   0.086562 112.273962  -0.365  0.71611    
## Light                0.005428   0.005991 113.810789   0.906  0.36688    
## DomD                 0.401155   0.054175 113.136002   7.405 2.52e-11 ***
## Seasonwinter:Light   0.029607   0.010208 113.291344   2.900  0.00448 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light  DomD  
## Seasonwintr -0.439                     
## Light       -0.540  0.593              
## DomD         0.099 -0.153 -0.393       
## Ssnwntr:Lgh  0.247 -0.817 -0.417 -0.001
ranova(add)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Y ~ Season + Light + Dom + (1 | Location) + Season:Light
##                npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>            7 -15.182 44.364                         
## (1 | Location)    6 -26.595 65.190 22.826  1  1.773e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(anova(add, type=2), digits=5)
## Type II Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Season       0.00760 0.00760     1 112.27  0.1329  0.716111    
## Light        0.39237 0.39237     1 114.92  6.8662  0.009971 ** 
## Dom          3.13331 3.13331     1 113.14 54.8302 2.521e-11 ***
## Season:Light 0.48067 0.48067     1 113.29  8.4113  0.004480 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.symb=var.table
var.symb$response<-"symb"

posthoc<-emmeans(add, ~Dom)
CLD(posthoc, Letters=letters)
##  Dom   emmean         SE   df lower.CL upper.CL .group
##  C   14.45561 0.08380064 3.26 14.20043 14.71080  a    
##  D   14.85677 0.09108456 4.53 14.61520 15.09834   b   
## 
## Results are averaged over the levels of: Season 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(add), ylab="symbiont density", par.strip.text=list(cex=0.7))

Chlorophyll a

######### chltotal-- 
Y<-model.data$chltot
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+Season:Dom + (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model
## Data: model.data
## Models:
## add: Y ~ Season + Light + Dom + Season:Dom + (1 | Location)
## full: Y ~ Season * Light * Dom + (1 | Location)
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## add   7 424.92 444.44 -205.46   410.92                         
## full 10 425.73 453.61 -202.87   405.73 5.1938      3     0.1581
summary(add)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Y ~ Season + Light + Dom + Season:Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 418.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.04486 -0.63880 -0.07154  0.47425  2.98104 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.2169   0.4657  
##  Residual             1.7844   1.3358  
## Number of obs: 120, groups:  Location, 4
## 
## Fixed effects:
##                    Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)         5.84785    0.40030  13.07577  14.609 1.77e-09 ***
## Seasonwinter        1.70590    0.30499 113.70630   5.593 1.56e-07 ***
## Light              -0.09232    0.03079 109.93020  -2.998  0.00336 ** 
## DomD               -0.49207    0.44865 114.16680  -1.097  0.27504    
## Seasonwinter:DomD  -1.32504    0.56779 113.45499  -2.334  0.02137 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light  DomD  
## Seasonwintr -0.528                     
## Light       -0.644  0.315              
## DomD         0.074  0.148 -0.459       
## Ssnwntr:DmD  0.017 -0.414  0.247 -0.741
ranova(add)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Y ~ Season + Light + Dom + (1 | Location) + Season:Dom
##                npar  logLik    AIC    LRT Df Pr(>Chisq)  
## <none>            7 -209.12 432.24                       
## (1 | Location)    6 -211.55 435.10 4.8651  1    0.02741 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(anova(add, type=2), digits=5)
## Type II Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Season     45.581  45.581     1 114.27 25.5446 1.658e-06 ***
## Light      16.042  16.042     1 109.93  8.9901  0.003357 ** 
## Dom        31.622  31.622     1 114.55 17.7213 5.111e-05 ***
## Season:Dom  9.718   9.718     1 113.45  5.4461  0.021372 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.chl=var.table
var.chl$response<-"chl"

posthoc<-emmeans(add, ~Season)
CLD(posthoc, Letters=letters)
##  Season   emmean        SE   df lower.CL upper.CL .group
##  summer 4.865228 0.3288109 6.70 4.080513 5.649944  a    
##  winter 5.908608 0.3026093 4.98 5.129885 6.687331   b   
## 
## Results are averaged over the levels of: Dom 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-emmeans(add, ~Dom|Season)
CLD(posthoc, Letters=letters)
## Season = summer:
##  Dom   emmean        SE    df lower.CL upper.CL .group
##  D   4.619195 0.4738171 24.45 3.642243 5.596147  a    
##  C   5.111262 0.3066745  5.30 4.336006 5.886518  a    
## 
## Season = winter:
##  Dom   emmean        SE    df lower.CL upper.CL .group
##  D   5.000055 0.3888373 12.97 4.159808 5.840303  a    
##  C   6.817161 0.3257516  6.50 6.034599 7.599723   b   
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(add), ylab="total chlor", par.strip.text=list(cex=0.7))

Chlorophyll per symbiont cell

######### chlcell -- 
Y<-model.data$chlcell
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model
## Data: model.data
## Models:
## add: Y ~ Season + Light + Dom + (1 | Location)
## full: Y ~ Season * Light * Dom + (1 | Location)
##      Df    AIC    BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## add   6 4.3051 21.030 3.8475  -7.6949                         
## full 10 6.3305 34.205 6.8347 -13.6695 5.9745      4     0.2011
summary(add)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 13.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.73631 -0.52024 -0.08737  0.53912  2.46089 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.02493  0.1579  
##  Residual             0.05191  0.2278  
## Number of obs: 120, groups:  Location, 4
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)    1.874743   0.096951   5.720514  19.337 1.98e-06 ***
## Seasonwinter   0.079964   0.047578 113.863005   1.681   0.0956 .  
## Light         -0.035217   0.005194 115.720118  -6.781 5.33e-10 ***
## DomD          -0.513165   0.051649 114.025011  -9.936  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light 
## Seasonwintr -0.411              
## Light       -0.480  0.481       
## DomD         0.100 -0.267 -0.433
ranova(add)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Y ~ Season + Light + Dom + (1 | Location)
##                npar   logLik    AIC    LRT Df Pr(>Chisq)    
## <none>            6  -6.5893 25.179                         
## (1 | Location)    5 -18.4197 46.839 23.661  1  1.149e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(anova(add, type=2), digits=5)
## Type II Analysis of Variance Table with Satterthwaite's method
##        Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Season 0.1466  0.1466     1 113.86  2.8247   0.09556 .  
## Light  2.3869  2.3869     1 115.72 45.9772 5.331e-10 ***
## Dom    5.1248  5.1248     1 114.03 98.7155 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.chlcell=var.table
var.chlcell$response<-"chlcell"

plot(allEffects(add), ylab="chla cell", par.strip.text=list(cex=0.7))

#ranef(add)
#ranova(add)
#fixef(add)
#sjp.lmer(add, y.offset = .4)

Isotope models

host d13C

########## host..d13C -- 
Y<-model.data$host..d13C
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ Season:Dom+ Season:Light+(1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model
## Data: model.data
## Models:
## add: Y ~ Season + Light + Dom + Season:Dom + Season:Light + (1 | Location)
## full: Y ~ Season * Light * Dom + (1 | Location)
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## add   8 323.78 346.08 -153.89   307.78                         
## full 10 327.70 355.58 -153.85   307.70 0.0775      2      0.962
summary(add)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## Y ~ Season + Light + Dom + Season:Dom + Season:Light + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 323.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.77180 -0.68022 -0.08845  0.75170  2.22371 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.2658   0.5156  
##  Residual             0.7386   0.8594  
## Number of obs: 120, groups:  Location, 4
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)        -16.41116    0.34694   7.99485 -47.303 4.47e-11 ***
## Seasonwinter         0.01219    0.31119 111.35688   0.039   0.9688    
## Light                0.13233    0.02315 112.82092   5.717 9.01e-08 ***
## DomD                -1.62857    0.29962 111.58775  -5.436 3.25e-07 ***
## Seasonwinter:DomD    0.84855    0.38854 111.42769   2.184   0.0311 *  
## Seasonwinter:Light  -0.01870    0.03892 112.30031  -0.480   0.6319    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light  DomD   Ssn:DD
## Seasonwintr -0.463                            
## Light       -0.558  0.557                     
## DomD         0.130 -0.110 -0.517              
## Ssnwntr:DmD -0.083  0.014  0.368 -0.760       
## Ssnwntr:Lgh  0.274 -0.775 -0.491  0.254 -0.335
ranova(add)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Y ~ Season + Light + Dom + (1 | Location) + Season:Dom + Season:Light
##                npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>            8 -161.60 339.20                         
## (1 | Location)    7 -171.63 357.27 20.063  1  7.494e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(anova(add), digits=5)
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Season        1.0324  1.0324     1 111.22  1.3978   0.23961    
## Light        23.6384 23.6384     1 113.94 32.0053 1.160e-07 ***
## Dom          27.4555 27.4555     1 112.24 37.1735 1.566e-08 ***
## Season:Dom    3.5227  3.5227     1 111.43  4.7696   0.03106 *  
## Season:Light  0.1704  0.1704     1 112.30  0.2307   0.63192    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.d13C.H=var.table
var.d13C.H$response<-"d13C.H"

posthoc<-emmeans(add, ~Dom)
CLD(posthoc, Letters=letters)
##  Dom    emmean        SE   df  lower.CL  upper.CL .group
##  D   -16.62814 0.3090539 5.04 -17.42049 -15.83578  a    
##  C   -15.42384 0.2779425 3.31 -16.26338 -14.58431   b   
## 
## Results are averaged over the levels of: Season 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-emmeans(add, ~Dom| Season)
CLD(posthoc, Letters=letters)
## Season = summer:
##  Dom    emmean        SE    df  lower.CL  upper.CL .group
##  D   -16.98392 0.3780843 10.98 -17.81631 -16.15154  a    
##  C   -15.35535 0.2880146  3.83 -16.16948 -14.54121   b   
## 
## Season = winter:
##  Dom    emmean        SE    df  lower.CL  upper.CL .group
##  D   -16.27236 0.3270155  6.30 -17.06346 -15.48125  a    
##  C   -15.49234 0.3058112  4.79 -16.28894 -14.69573   b   
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(add), ylab="d13Chost", par.strip.text=list(cex=0.7))

symbiont d13C

########## symb..d13C --
Y<-model.data$symb..d13C
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ Season:Dom +(1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model
## Data: model.data
## Models:
## add: Y ~ Season + Light + Dom + Season:Dom + (1 | Location)
## full: Y ~ Season * Light * Dom + (1 | Location)
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## add   7 334.15 353.66 -160.08   320.15                         
## full 10 339.80 367.68 -159.90   319.80 0.3464      3     0.9511
summary(add)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Y ~ Season + Light + Dom + Season:Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 330.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9724 -0.6051 -0.0345  0.7590  2.8420 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.3887   0.6234  
##  Residual             0.8043   0.8969  
## Number of obs: 120, groups:  Location, 4
## 
## Fixed effects:
##                    Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)       -16.59013    0.38242   5.80200 -43.382 1.64e-08 ***
## Seasonwinter       -0.29714    0.20544 112.63287  -1.446   0.1509    
## Light               0.14086    0.02111 114.73832   6.673 9.30e-10 ***
## DomD               -1.53545    0.30257 112.84432  -5.075 1.54e-06 ***
## Seasonwinter:DomD   1.26115    0.38222 112.51114   3.300   0.0013 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light  DomD  
## Seasonwintr -0.377                     
## Light       -0.462  0.322              
## DomD         0.061  0.141 -0.466       
## Ssnwntr:DmD  0.009 -0.411  0.249 -0.741
ranova(add)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Y ~ Season + Light + Dom + (1 | Location) + Season:Dom
##                npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>            7 -165.09 344.17                         
## (1 | Location)    6 -178.53 369.07 26.894  1  2.149e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(anova(add, type=2), digits=5)
## Type II Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Season      0.002   0.002     1 112.90  0.0023 0.9615036    
## Light      35.816  35.816     1 114.74 44.5288 9.298e-10 ***
## Dom        12.375  12.375     1 113.05 15.3856 0.0001508 ***
## Season:Dom  8.757   8.757     1 112.51 10.8868 0.0012975 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.d13C.S=var.table
var.d13C.S$response<-"d13C.S"

posthoc<-emmeans(add, ~Dom)
CLD(posthoc, Letters=letters)
##  Dom    emmean        SE   df  lower.CL  upper.CL .group
##  D   -16.51965 0.3551755 4.34 -17.47613 -15.56317  a    
##  C   -15.61478 0.3274383 3.15 -16.62935 -14.60020   b   
## 
## Results are averaged over the levels of: Season 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-emmeans(add, ~Dom| Season)
CLD(posthoc, Letters=letters)
## Season = summer:
##  Dom    emmean        SE   df  lower.CL  upper.CL .group
##  D   -17.00166 0.4167784 8.07 -17.96137 -16.04195  a    
##  C   -15.46621 0.3392889 3.63 -16.44701 -14.48541   b   
## 
## Season = winter:
##  Dom    emmean        SE   df  lower.CL  upper.CL .group
##  D   -16.03764 0.3752434 5.40 -16.98092 -15.09436  a    
##  C   -15.76334 0.3471338 3.96 -16.73061 -14.79608  a    
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(add), ylab="d13Csymb", par.strip.text=list(cex=0.7))

host-symbiont d13C

######### d13C..host.sym --
Y<-model.data$d13C..host.sym
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ Season:Light +Season:Dom +(1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model
## Data: model.data
## Models:
## add: Y ~ Season + Light + Dom + Season:Light + Season:Dom + (1 | Location)
## full: Y ~ Season * Light * Dom + (1 | Location)
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## add   8 116.04 138.34 -50.022  100.044                         
## full 10 119.77 147.64 -49.884   99.768 0.2757      2     0.8712
summary(add)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## Y ~ Season + Light + Dom + Season:Light + Season:Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 126.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.9375 -0.3451  0.0104  0.4548  3.6783 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.02824  0.1680  
##  Residual             0.13287  0.3645  
## Number of obs: 120, groups:  Location, 4
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)   
## (Intercept)          0.126970   0.129306  10.967796   0.982  0.34730   
## Seasonwinter         0.409353   0.131934 111.362679   3.103  0.00243 **
## Light               -0.002278   0.009790 113.441986  -0.233  0.81639   
## DomD                -0.136399   0.126986 111.723004  -1.074  0.28508   
## Seasonwinter:Light  -0.034256   0.016478 112.783500  -2.079  0.03990 * 
## Seasonwinter:DomD   -0.347111   0.164712 111.486057  -2.107  0.03733 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light  DomD   Ssnw:L
## Seasonwintr -0.528                            
## Light       -0.633  0.560                     
## DomD         0.147 -0.111 -0.516              
## Ssnwntr:Lgh  0.315 -0.775 -0.497  0.257       
## Ssnwntr:DmD -0.094  0.015  0.370 -0.761 -0.336
ranova(add)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Y ~ Season + Light + Dom + (1 | Location) + Season:Light + Season:Dom
##                npar  logLik    AIC    LRT Df Pr(>Chisq)   
## <none>            8 -63.134 142.27                        
## (1 | Location)    7 -67.238 148.48 8.2073  1   0.004172 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(anova(add, type=2), digits=4)
## Type II Analysis of Variance Table with Satterthwaite's method
##              Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)    
## Season       1.3196  1.3196     1 111.4   9.931  0.00209 ** 
## Light        0.3604  0.3604     1 112.9   2.712  0.10235    
## Dom          2.2911  2.2911     1 112.8  17.243 6.42e-05 ***
## Season:Light 0.5742  0.5742     1 112.8   4.322  0.03990 *  
## Season:Dom   0.5901  0.5901     1 111.5   4.441  0.03733 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.d13C.HS=var.table
var.d13C.HS$response<-"d13C.HS"

posthoc<-emmeans(add, ~Dom)
CLD(posthoc, Letters=letters)
##  Dom     emmean         SE   df   lower.CL  upper.CL .group
##  D   -0.1331473 0.11086251 6.46 -0.3998082 0.1335136  a    
##  C    0.1768072 0.09489086 3.49 -0.1025094 0.4561238   b   
## 
## Results are averaged over the levels of: Season 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-emmeans(add, ~Dom|Season)
CLD(posthoc, Letters=letters)
## Season = summer:
##  Dom      emmean        SE    df   lower.CL  upper.CL .group
##  D   -0.02760769 0.1443214 17.32 -0.3316682 0.2764528  a    
##  C    0.10879151 0.1001424  4.36 -0.1604178 0.3780008  a    
## 
## Season = winter:
##  Dom      emmean        SE    df   lower.CL  upper.CL .group
##  D   -0.23868691 0.1197827  8.71 -0.5110137 0.0336399  a    
##  C    0.24482291 0.1092528  5.96 -0.0229421 0.5125880   b   
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(add), ylab="d13C-hs", par.strip.text=list(cex=0.7))

skeleton d13C

####### d13C..skel --
Y<-model.data$d13C..skel
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model
## Data: model.data
## Models:
## add: Y ~ Season + Light + Dom + (1 | Location)
## full: Y ~ Season * Light * Dom + (1 | Location)
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## add   6 313.69 330.42 -150.85   301.69                         
## full 10 319.03 346.91 -149.52   299.03 2.6592      4     0.6164
summary(add)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 313
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.96117 -0.81100 -0.02084  0.65228  2.40914 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.1400   0.3742  
##  Residual             0.7022   0.8380  
## Number of obs: 120, groups:  Location, 4
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   -2.855001   0.277939   9.891174 -10.272 1.36e-06 ***
## Seasonwinter   0.460523   0.174545 114.750779   2.638  0.00949 ** 
## Light          0.008887   0.018898 115.177685   0.470  0.63906    
## DomD           0.011169   0.189383 115.013233   0.059  0.95308    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light 
## Seasonwintr -0.522              
## Light       -0.610  0.477       
## DomD         0.122 -0.263 -0.428
ranova(add)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Y ~ Season + Light + Dom + (1 | Location)
##                npar  logLik    AIC    LRT Df Pr(>Chisq)   
## <none>            6 -156.48 324.97                        
## (1 | Location)    5 -161.74 333.47 10.506  1    0.00119 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(anova(add, type=2), digits=4)
## Type II Analysis of Variance Table with Satterthwaite's method
##        Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)   
## Season  4.888   4.888     1 114.8   6.961 0.00949 **
## Light   0.155   0.155     1 115.2   0.221 0.63906   
## Dom     0.002   0.002     1 115.0   0.003 0.95308   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.d13C.sk=var.table
var.d13C.sk$response<-"d13C.sk"

posthoc<-emmeans(add, ~Season)
CLD(posthoc, Letters=letters)
##  Season    emmean        SE   df  lower.CL  upper.CL .group
##  summer -2.778508 0.2286013 4.69 -3.377847 -2.179170  a    
##  winter -2.317985 0.2197840 4.07 -2.924134 -1.711836   b   
## 
## Results are averaged over the levels of: Dom 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(add), ylab="d13C-skel", par.strip.text=list(cex=0.7))

host d15N

####### host..d15N --
Y<-model.data$host..d15N
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+(1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model
## Data: model.data
## Models:
## add: Y ~ Season + Light + Dom + (1 | Location)
## full: Y ~ Season * Light * Dom + (1 | Location)
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## add   6 86.080 102.81 -37.040   74.080                         
## full 10 88.549 116.42 -34.275   68.549 5.5304      4     0.2371
summary(add)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 90.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4078 -0.5566  0.0705  0.6116  2.4030 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.29199  0.5404  
##  Residual             0.09669  0.3110  
## Number of obs: 120, groups:  Location, 4
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)    4.867242   0.280992   3.425434  17.322 0.000186 ***
## Seasonwinter   0.069196   0.065050 113.164781   1.064 0.289709    
## Light         -0.014857   0.007142 113.658553  -2.080 0.039755 *  
## DomD           0.044236   0.070641 113.196111   0.626 0.532435    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light 
## Seasonwintr -0.195              
## Light       -0.228  0.484       
## DomD         0.049 -0.270 -0.436
ranova(add)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Y ~ Season + Light + Dom + (1 | Location)
##                npar   logLik    AIC LRT Df Pr(>Chisq)    
## <none>            6  -45.329 102.66                      
## (1 | Location)    5 -105.829 221.66 121  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(anova(add, type=2), digits=4)
## Type II Analysis of Variance Table with Satterthwaite's method
##        Sum Sq Mean Sq NumDF DenDF F value Pr(>F)  
## Season 0.1094  0.1094     1 113.2   1.132 0.2897  
## Light  0.4184  0.4184     1 113.7   4.327 0.0398 *
## Dom    0.0379  0.0379     1 113.2   0.392 0.5324  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.d15N.H=var.table
var.d15N.H$response<-"d15N.H"

plot(allEffects(add), ylab="d15Nhost", par.strip.text=list(cex=0.7))

symbiont d15N

######## symb..d15N --
Y<-model.data$symb..d15N
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ Season:Dom+ Light:Dom + (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model
## Data: model.data
## Models:
## add: Y ~ Season + Light + Dom + Season:Dom + Light:Dom + (1 | Location)
## full: Y ~ Season * Light * Dom + (1 | Location)
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## add   8 103.85 126.15 -43.925   87.850                         
## full 10 107.53 135.40 -43.763   87.526 0.3242      2     0.8504
summary(add)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## Y ~ Season + Light + Dom + Season:Dom + Light:Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 112.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.10105 -0.57497  0.06009  0.58736  2.67836 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.4296   0.6555  
##  Residual             0.1094   0.3307  
## Number of obs: 120, groups:  Location, 4
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)         4.343041   0.339414   3.386078  12.796 0.000551 ***
## Seasonwinter       -0.080336   0.076829 111.087667  -1.046 0.297994    
## Light              -0.004405   0.008747 111.423764  -0.504 0.615491    
## DomD               -0.183898   0.247719 111.041211  -0.742 0.459434    
## Seasonwinter:DomD   0.435748   0.179700 111.065898   2.425 0.016926 *  
## Light:DomD          0.013734   0.016369 111.016288   0.839 0.403249    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light  DomD   Ssn:DD
## Seasonwintr -0.171                            
## Light       -0.216  0.358                     
## DomD        -0.073  0.203  0.204              
## Ssnwntr:DmD  0.061 -0.416 -0.098 -0.815       
## Light:DomD   0.095 -0.158 -0.442 -0.892  0.619
ranova(add)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Y ~ Season + Light + Dom + (1 | Location) + Season:Dom + Light:Dom
##                npar   logLik    AIC    LRT Df Pr(>Chisq)    
## <none>            8  -56.105 128.21                         
## (1 | Location)    7 -125.714 265.43 139.22  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(anova(add, type=2), digits=6)
## Type II Analysis of Variance Table with Satterthwaite's method
##              Sum Sq  Mean Sq NumDF   DenDF F value   Pr(>F)  
## Season     0.000027 0.000027     1 111.116 0.00025 0.987387  
## Light      0.002411 0.002411     1 111.517 0.02204 0.882246  
## Dom        0.506556 0.506556     1 111.038 4.63087 0.033568 *
## Season:Dom 0.643192 0.643192     1 111.066 5.87997 0.016926 *
## Light:Dom  0.077007 0.077007     1 111.016 0.70399 0.403249  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.d15N.S=var.table
var.d15N.S$response<-"d15N.S"

posthoc<-emmeans(add, ~Dom|Season)
CLD(posthoc, Letters=letters)
## Season = summer:
##  Dom   emmean        SE   df lower.CL upper.CL .group
##  D   4.233574 0.3545527 4.03 3.251927 5.215221  a    
##  C   4.307890 0.3314409 3.08 3.268355 5.347426  a    
## 
## Season = winter:
##  Dom   emmean        SE   df lower.CL upper.CL .group
##  C   4.227554 0.3327041 3.13 3.192580 5.262528  a    
##  D   4.588986 0.3368169 3.28 3.567657 5.610315   b   
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(add), ylab="d15Nsymb", par.strip.text=list(cex=0.7))

host-symbiont d15N

####### d15N..host.sym  -- 
Y<-model.data$d15N..host.sym
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ Season:Dom+  (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model
## Data: model.data
## Models:
## add: Y ~ Season + Light + Dom + Season:Dom + (1 | Location)
## full: Y ~ Season * Light * Dom + (1 | Location)
##      Df    AIC    BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## add   7 11.866 31.378 1.0670  -2.1341                         
## full 10 15.598 43.472 2.2013  -4.4025 2.2685      3     0.5186
summary(add)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Y ~ Season + Light + Dom + Season:Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 21.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0812 -0.6452  0.1012  0.7207  2.0871 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.01312  0.1145  
##  Residual             0.05606  0.2368  
## Number of obs: 120, groups:  Location, 4
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)         0.538132   0.081629   8.560016   6.592 0.000126 ***
## Seasonwinter        0.162519   0.054163 113.073851   3.001 0.003316 ** 
## Light              -0.013270   0.005526 114.609281  -2.401 0.017938 *  
## DomD                0.078624   0.079731 113.429934   0.986 0.326175    
## Seasonwinter:DomD  -0.417705   0.100797 112.873708  -4.144 6.62e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light  DomD  
## Seasonwintr -0.463                     
## Light       -0.566  0.319              
## DomD         0.071  0.144 -0.463       
## Ssnwntr:DmD  0.012 -0.412  0.248 -0.741
ranova(add)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Y ~ Season + Light + Dom + (1 | Location) + Season:Dom
##                npar  logLik    AIC    LRT Df Pr(>Chisq)   
## <none>            7 -10.947 35.895                        
## (1 | Location)    6 -16.249 44.498 10.603  1   0.001129 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(anova(add, type=2), digits=4)
## Type II Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)    
## Season     0.1037  0.1037     1 113.5   1.849  0.17659    
## Light      0.3233  0.3233     1 114.6   5.767  0.01794 *  
## Dom        0.5375  0.5375     1 113.8   9.588  0.00247 ** 
## Season:Dom 0.9627  0.9627     1 112.9  17.173 6.62e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.d15N.HS<-var.table
var.d15N.HS$response<-"d15N.HS"

posthoc<-emmeans(add, ~Dom|Season)
CLD(posthoc, Letters=letters)
## Season = summer:
##  Dom    emmean         SE    df  lower.CL  upper.CL .group
##  C   0.4322534 0.06730944  4.27 0.2498890 0.6146178  a    
##  D   0.5108774 0.09284886 14.23 0.3120315 0.7097232  a    
## 
## Season = winter:
##  Dom    emmean         SE    df  lower.CL  upper.CL .group
##  D   0.2556905 0.07951643  8.15 0.0729050 0.4384760  a    
##  C   0.5947720 0.07005404  4.94 0.4139888 0.7755552   b   
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(add), ylab="d15N-hs", par.strip.text=list(cex=0.7))

host C:N

###### host..C.N --
Y<-model.data$host..C.N
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model
## Data: model.data
## Models:
## add: Y ~ Season + Light + Dom + (1 | Location)
## full: Y ~ Season * Light * Dom + (1 | Location)
##      Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## add   6 -273.59 -256.87 142.79  -285.59                         
## full 10 -268.55 -240.68 144.28  -288.55 2.9622      4     0.5642
summary(add)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: -254.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7879 -0.7136 -0.1026  0.6668  3.2204 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Location (Intercept) 0.0008847 0.02974 
##  Residual             0.0052870 0.07271 
## Number of obs: 120, groups:  Location, 4
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  1.944e+00  2.319e-02 1.050e+01  83.819 3.46e-16 ***
## Seasonwinter 3.999e-03  1.513e-02 1.149e+02   0.264    0.792    
## Light        1.371e-03  1.635e-03 1.142e+02   0.839    0.403    
## DomD         6.456e-03  1.642e-02 1.152e+02   0.393    0.695    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light 
## Seasonwintr -0.541              
## Light       -0.632  0.476       
## DomD         0.125 -0.262 -0.427
ranova(add)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Y ~ Season + Light + Dom + (1 | Location)
##                npar logLik     AIC    LRT Df Pr(>Chisq)   
## <none>            6 127.30 -242.59                        
## (1 | Location)    5 123.86 -237.72 6.8693  1   0.008769 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(anova(add, type=2), digits=4)
## Type II Analysis of Variance Table with Satterthwaite's method
##          Sum Sq  Mean Sq NumDF DenDF F value Pr(>F)
## Season 0.000369 0.000369     1 114.9   0.070  0.792
## Light  0.003718 0.003718     1 114.2   0.703  0.403
## Dom    0.000818 0.000818     1 115.2   0.155  0.695
var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.CN.H<-var.table
var.CN.H$response<-"CN.Host"

plot(allEffects(add), ylab="C:N host", par.strip.text=list(cex=0.7))

symbiont C:N

####### symb..C.N -- 
Y<-model.data$symb..C.N
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model
## Data: model.data
## Models:
## add: Y ~ Season + Light + Dom + (1 | Location)
## full: Y ~ Season * Light * Dom + (1 | Location)
##      Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## add   6 -212.29 -195.61 112.14  -224.29                         
## full 10 -207.16 -179.37 113.58  -227.16 2.8707      4     0.5797
summary(add)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: -193.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8398 -0.6571 -0.1712  0.6883  2.9658 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.000000 0.00000 
##  Residual             0.009201 0.09592 
## Number of obs: 119, groups:  Location, 4
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   2.051e+00  2.207e-02  1.150e+02  92.903   <2e-16 ***
## Seasonwinter -2.956e-02  1.957e-02  1.150e+02  -1.510    0.134    
## Light         3.778e-04  1.953e-03  1.150e+02   0.193    0.847    
## DomD          8.659e-04  2.110e-02  1.150e+02   0.041    0.967    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light 
## Seasonwintr -0.684              
## Light       -0.798  0.433       
## DomD         0.097 -0.230 -0.379
## convergence code: 0
## singular fit
ranova(add)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Y ~ Season + Light + Dom + (1 | Location)
##                npar logLik     AIC         LRT Df Pr(>Chisq)
## <none>            6 96.835 -181.67                          
## (1 | Location)    5 96.835 -183.67 -2.8422e-14  1          1
print(anova(add, type=2), digits=5)
## Type II Analysis of Variance Table with Satterthwaite's method
##           Sum Sq   Mean Sq NumDF DenDF F value Pr(>F)
## Season 0.0209902 0.0209902     1   115  2.2813 0.1337
## Light  0.0003443 0.0003443     1   115  0.0374 0.8470
## Dom    0.0000155 0.0000155     1   115  0.0017 0.9673
var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.CN.S<-var.table
var.CN.S$response<-"CN.Symb"

plot(allEffects(add), ylab="C:N symb", par.strip.text=list(cex=0.7))

Using compiled models above, a summary table and figure for random effects can be generated. This table figure shows the proportion of variance explained by the random effect Location. For symbiont tissue C:N, Location explained <0.00001 of variance.

# random effect table

rand.eff.table<-rbind(var.biomass, var.symb, var.chl, var.chlcell, var.d13C.H, var.d13C.S, var.d13C.HS, var.d13C.sk, var.d15N.H, var.d15N.S, var.d15N.HS, var.CN.H, var.CN.S); rand.eff.table=rand.eff.table[-2:-3]
colnames(rand.eff.table)<-c("group", "variance", "stdev", "prop.var", "response")
rand.eff.table=rand.eff.table[,c(5,1:4)]
rand.eff.summary<-rand.eff.table[,c(1,5)]

rand.eff.df<- rand.eff.summary[-which(rand.eff.summary$prop.var==""), ] # removing blank values
rand.eff.df$prop.var<-as.numeric(rand.eff.df$prop.var)

rand.eff.df$response<-as.factor(rand.eff.df$response)
rand.eff.df$response<-factor(rand.eff.df$response, levels = c("biomass", "symb", "chl", "chlcell", 
                                        "d13C.H", "d13C.S", "d13C.HS", "d13C.sk", 
                                        "d15N.H", "d15N.S", "d15N.HS",
                                        "CN.Host", "CN.Symb"))

var.labs<-c("biomass", "symbionts", "chlorophyll", "chl/cell", 
            (expression(paste(delta^{13}, C[H]))), (expression(paste(delta^{13}, C[S]))),
            (expression(paste(delta^{13}, C[H-S]))), (expression(paste(delta^{13}, C[Sk]))),
            (expression(paste(delta^{15}, N[H]))), (expression(paste(delta^{15}, N[S]))),
            (expression(paste(delta^{15}, N[H-S]))), 
            (expression(paste(C:N[H]))), (expression(paste(C:N[S]))))


# make plot
Fig.rand.eff<-ggplot(rand.eff.df, aes(x=response, y=prop.var, fill=response)) +
  geom_bar(stat="identity", position = pd, width=0.7, size=0.4) +
  theme(axis.text.x = element_text(size=10, angle=45, hjust = 1)) +
  scale_x_discrete(name="Responses", label=var.labs) +
  theme(axis.title.y= element_text(size=10),
        axis.text.y= element_text(size=10))+
  ylab(expression(paste("variance proportion", sep=""))) +
  scale_y_continuous(expand=c(0,0), limits=c(0, 1)) + theme(legend.position="none"); Fig.rand.eff
**Figure** *Proportion of model variance explained by the random effect of **Location**, representing the four reef locations where corals were collected within each season*

Figure Proportion of model variance explained by the random effect of Location, representing the four reef locations where corals were collected within each season 

ggsave(file="figures/models/variance.pdf", height=3.5, width=6)

Figures

Physiology

#########################
# all data
df.fig<-data.trim[,c(17, 1:5, 18:25, 7:16, 26:29)]
table(df.fig$Dom:df.fig$Season)

# figure layout
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

##################
################## 

## season relationship
data.summer<-df.fig[df.fig$Season=="summer", ]
data.winter<-df.fig[df.fig$Season=="winter", ]

## Symbiont relationship
C.sum.df<-df.fig[(df.fig$Dom=="C" & df.fig$Season=="summer") ,]
C.win.df<-df.fig[(df.fig$Dom=="C" & df.fig$Season=="winter") ,]

D.sum.df<-df.fig[(df.fig$Dom=="D" & df.fig$Season=="summer") ,]
D.win.df<-df.fig[(df.fig$Dom=="D" & df.fig$Season=="winter") ,]

##transformations
# model.data$zoox<-log(model.data$zoox)
# model.data$chlcell<-log(model.data$chlcell)
# model.data$host..C.N<-log(model.data$host..C.N)
# model.data$symb..C.N<-log(model.data$symb..C.N)

Tissue biomass

##################
# Fig: biomass over season and depth
################## 
par(mfrow=c(1,2),mar=c(5,4.5,2,2))
mod<-lmer(biomass~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy_lmer(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params

# C summer
plot(biomass~Light, data=C.sum.df, col='red3', ylim=c(0, 80), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon", x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
par(new=T)

# D summer
plot(biomass~Light, data=D.sum.df, col='blue', ylim=c(0, 80), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab=expression(paste("total chlorophyll", ~(mu*g~cm^-2), sep="")), cex.axis=0.8,
     cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]), coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)

# C winter
plot(biomass~Light, data=C.win.df, col='red3', ylim=c(0, 80), xlim=c(15, 0), yaxt="n", xaxt="n",
     cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]), coef(summary(mod))["Light",1], col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)

# D winter
plot(biomass~Light, data=D.win.df, col='blue', ylim=c(0, 80), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab="", cex.axis=0.8, cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1]) , coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)

dev.copy(pdf, "figures/physiology/PanKB.biomass.pdf", width=6, height=4)
dev.off()

Total chlorophyll

##################
# Fig: chlorophyll (total) over season and depth
################## 
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# model 
mod<-lmer(chltot~Season+Light+Dom+Season:Dom + (1|Location), data=model.data, na.action=na.exclude)
fixed_params <- tidy_lmer(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params

intercept<-coef(summary(mod))[1] # intercept
coef(summary(mod))["Seasonwinter",1] # winter, summer = 0
coef(summary(mod))["Light",1] # light relationship (slope for all figures)
coef(summary(mod))["DomD",1] # domD, domC = 0
coef(summary(mod))["Seasonwinter:DomD",1] #D in winter, winter C = 0, summer C = 0

# C summer
plot(chltot~Light, data=C.sum.df, col='red3', ylim=c(0, 15), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon", x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
par(new=T)

# D summer
plot(chltot~Light, data=D.sum.df, col='blue', ylim=c(0, 15), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab=expression(paste("total chlorophyll", ~(mu*g~cm^-2), sep="")), cex.axis=0.8,
     cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]), coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)

# C winter
plot(chltot~Light, data=C.win.df, col='red3', ylim=c(0, 15), xlim=c(15, 0), yaxt="n", xaxt="n",cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]), coef(summary(mod))["Light",1], col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)

# D winter
plot(chltot~Light, data=D.win.df, col='blue', ylim=c(0, 15), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab="",
     cex.lab=0.8, cex.axis=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1] + coef(summary(mod))["Seasonwinter:DomD",1]) , coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)

dev.copy(pdf, "figures/physiology/PanKB.chla.pdf", width=6, height=4)
dev.off()

Chlorophyll per symbiont cell

##################
# Fig: chlorophyll/cell over season and depth
##################  
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# model.data$chlcell<-log(model.data$chlacell)
mod<-lmer(chlcell~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy_lmer(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params

##### back transformed
x_light=seq(1,23, 0.5)
ypred= exp(intercept+ coef(summary(mod))["Light",1]*x_light)

# C summer
plot(chlcell~Light, data=C.sum.df, col='red3', ylim=c(0, 15), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8,
     pch=21, bg="salmon", ylab="", xlab ="")
lines(x_light, ypred, col="salmon", lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
par(new=T)

# D summer
x_light=seq(5,24, 0.5)
ypred= exp(intercept+ coef(summary(mod))["DomD",1]+coef(summary(mod))["Light",1]*x_light)

plot(chlcell~Light, data=D.sum.df, col='blue', ylim=c(0, 15), xlim=c(25, 0), cex=0.8, pch=21,
     bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
     ylab=(expression(paste("chlorophyll (",pg~cell^-1, ")", sep=""))), 
     cex.axis=0.8, cex.lab=0.8, main="Summer", cex.main=0.7)
lines(x_light, ypred, col="deepskyblue1", lwd=2)


# C winter
x_light=seq(0.3,14, 0.5)
ypred= exp(intercept + coef(summary(mod))["Seasonwinter",1] + coef(summary(mod))["Light",1]*x_light)

plot(chlcell~Light, data=C.win.df, col='red3', ylim=c(0, 15), xlim=c(15, 0), yaxt="n", xaxt="n",cex=0.8,
     pch=21, bg="salmon", ylab="", xlab ="")
lines(x_light, ypred, col="salmon", lwd=2)

par(new=T)

# D winter
x_light=seq(1.5,14, 0.5)
ypred= exp(intercept+ coef(summary(mod))["Seasonwinter",1] + coef(summary(mod))["DomD",1] +
             coef(summary(mod))["Light",1]*x_light)

plot(chlcell~Light, data=D.win.df, col='blue', ylim=c(0, 15), xlim=c(15, 0), cex=0.8, pch=21,
     bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
     ylab="",
     cex.axis=0.8, cex.lab=0.8, main="Winter", cex.main=0.7)
lines(x_light, ypred, col="deepskyblue1", lwd=2)

dev.copy(pdf, "figures/physiology/PanKB.chlacell.pdf", width=6, height=4)
dev.off()

Symbiont densities

##################
# Fig: symbionts over season and depth
##################  
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# model.data$zoox<-log(model.data$zoox)
mod<-lmer(zoox~Season+Light+Dom+Season:Light + (1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy_lmer(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params


####################    
##### back transformed
x_light=seq(1,23, 0.5)
ypred= exp(intercept+ coef(summary(mod))["Light",1]*x_light)

# C summer
plot(zoox~Light, data=C.sum.df, col='red3', ylim=c(0, 6000000), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="")
lines(x_light, ypred, col="salmon", lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
par(new=T)

# D summer
x_light=seq(5,24, 0.5)
ypred= exp(intercept+ coef(summary(mod))["DomD",1]+coef(summary(mod))["Light",1]*x_light)

plot(zoox~Light, data=D.sum.df, col='blue', ylim=c(0, 6000000), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab=expression(paste("symbionts ", "(",cells~cm^-2~")", sep="")), 
     cex.axis=0.8, cex.lab=0.8, main="Summer", cex.main=0.7)
lines(x_light, ypred, col="deepskyblue1", lwd=2)


# C winter
x_light=seq(0.3,14, 0.5)
ypred= exp(intercept + coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["Light",1]+ coef(summary(mod))["Seasonwinter:Light",1]*x_light)

plot(zoox~Light, data=C.win.df, col='red3', ylim=c(0, 6000000), xlim=c(15, 0), yaxt="n", xaxt="n",cex=0.8, pch=21, bg="salmon", ylab="", xlab ="")
lines(x_light, ypred, col="salmon", lwd=2)

par(new=T)

# D winter
x_light=seq(1.5,14, 0.5)
ypred= exp(intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["Light",1]+ coef(summary(mod))["DomD",1] + coef(summary(mod))["Seasonwinter:Light",1]*x_light)

plot(zoox~Light, data=D.win.df, col='blue', ylim=c(0, 6000000), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
     ylab="", cex.axis=0.8,
     cex.axis=0.8, cex.lab=0.8, main="Winter", cex.main=0.7)
lines(x_light, ypred, col="deepskyblue1", lwd=2)

dev.copy(pdf, "figures/physiology/PanKB.zoox.pdf", width=6, height=4)
dev.off()

Isotopes

Host tissue carbon isotope composition

##################
# Fig: d13C host over season and depth
##################        
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# model.data$host..d13C
mod<-lmer(host..d13C~Season+Light+Dom+ Season:Dom+ Season:Light+(1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy_lmer(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params

# C summer
plot(host..d13C~Light, data=C.sum.df, col='red3', ylim=c(-20, -10), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon", 
           x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
par(new=T)

# D summer
plot(host..d13C~Light, data=D.sum.df, col='blue', ylim=c(-20, -10), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab=(expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]), coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)

# C winter
plot(host..d13C~Light, data=C.win.df, col='red3', ylim=c(-20, -10), xlim=c(15, 0), yaxt="n", xaxt="n",
     cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]),   (coef(summary(mod))["Light",1]+coef(summary(mod))["Seasonwinter:Light",1]), col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)

# D winter
plot(host..d13C~Light, data=D.win.df, col='blue', ylim=c(-20, -10), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab="", cex.axis=0.8, cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1]+ coef(summary(mod))["Seasonwinter:DomD",1]), (coef(summary(mod))["Light",1]+coef(summary(mod))["Seasonwinter:Light",1]), 
      col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)

dev.copy(pdf, "figures/isotope/PanKB.d13C-host.pdf", width=6, height=4, encod="MacRoman")
dev.off()

Symbiont tissue carbon isotope composition

##################
# Fig: d13C symb over season and depth
##################        
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# model.data$symb..d13C
mod<-lmer(symb..d13C~Season+Light+Dom+ Season:Dom+ Season:Light+(1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy_lmer(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params

# C summer
plot(symb..d13C~Light, data=C.sum.df, col='red3', ylim=c(-20, -10), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon", 
           x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
par(new=T)

# D summer
plot(symb..d13C~Light, data=D.sum.df, col='blue', ylim=c(-20, -10), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab=(expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]), coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)

# C winter
plot(symb..d13C~Light, data=C.win.df, col='red3', ylim=c(-20, -10), xlim=c(15, 0), yaxt="n", xaxt="n",
     cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]), (coef(summary(mod))["Light",1]+ coef(summary(mod))["Seasonwinter:Light",1]), col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)

# D winter
plot(symb..d13C~Light, data=D.win.df, col='blue', ylim=c(-20, -10), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab="", cex.axis=0.8, cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1]+ coef(summary(mod))["Seasonwinter:DomD",1]),
(coef(summary(mod))["Light",1]+coef(summary(mod))["Seasonwinter:Light",1]), col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)

dev.copy(pdf, "figures/isotope/PanKB.d13C-symb.pdf", width=6, height=4, encod="MacRoman")
dev.off()

Host-symbiont tissue carbon isotope composition

##################
# Fig: d13C host-symb over season and depth
##################        
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# model.data$d13C..host.sym
mod<-lmer(d13C..host.sym~Season+Light+Dom+ Season:Light +Season:Dom +(1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy_lmer(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params

# C summer
plot(d13C..host.sym~Light, data=C.sum.df, col='red3', ylim=c(-3, 3), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon", 
           x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
par(new=T)

# D summer
plot(d13C..host.sym~Light, data=D.sum.df, col='blue', ylim=c(-3, 3), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab=(expression(paste(delta^{13}, C[H-S], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]), coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)

# C winter
plot(d13C..host.sym~Light, data=C.win.df, col='red3', ylim=c(-3, 3), xlim=c(15, 0), yaxt="n", xaxt="n",
     cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]), (coef(summary(mod))["Light",1]+coef(summary(mod))["Seasonwinter:Light",1]), col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)

# D winter
plot(d13C..host.sym~Light, data=D.win.df, col='blue', ylim=c(-3, 3), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab="", cex.axis=0.8, cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1]+
              coef(summary(mod))["Seasonwinter:DomD",1]),
           (coef(summary(mod))["Light",1]+ coef(summary(mod))["Seasonwinter:Light",1]), 
      col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)

dev.copy(pdf, "figures/isotope/PanKB.d13Ch-s.pdf", width=6, height=4, encod="MacRoman")
dev.off()

Coral skeleton tissue carbon isotope composition

##################
# Fig: d13C skeleton over season and depth
##################        
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# model.data$d13C..skel
mod<-lmer(d13C..skel~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy_lmer(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params

# C summer
plot(d13C..skel~Light, data=C.sum.df, col='red3', ylim=c(-8, 4), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon", 
           x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
par(new=T)

# D summer
plot(d13C..skel~Light, data=D.sum.df, col='blue', ylim=c(-8, 4), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab=(expression(paste(delta^{13}, C[Sk], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]), coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)

# C winter
plot(d13C..skel~Light, data=C.win.df, col='red3', ylim=c(-8, 4), xlim=c(15, 0), yaxt="n", xaxt="n",
     cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]), coef(summary(mod))["Light",1], col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)

# D winter
plot(d13C..skel~Light, data=D.win.df, col='blue', ylim=c(-8, 4), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab="", cex.axis=0.8, cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1]),
       coef(summary(mod))["Light",1], 
      col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)

dev.copy(pdf, "figures/isotope/PanKB.d13C-skel.pdf", width=6, height=4, encod="MacRoman")
dev.off()

Host tissue nitrogen isotope composition

##################
# Fig: d15N host over season and depth
##################        
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# model.data$host..d15N
mod<-lmer(host..d15N~Season+Light+Dom+(1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy_lmer(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params

# C summer
plot(host..d15N~Light, data=C.sum.df, col='red3', ylim=c(2, 8), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon", 
           x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
par(new=T)

# D summer
plot(host..d15N~Light, data=D.sum.df, col='blue', ylim=c(2, 8), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab=(expression(paste(delta^{15}, N[H], " (\u2030, Air)"))), cex.axis=0.8,
     cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]), coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)

# C winter
plot(host..d15N~Light, data=C.win.df, col='red3', ylim=c(2, 8), xlim=c(15, 0), yaxt="n", xaxt="n",
     cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]), coef(summary(mod))["Light",1], col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)

# D winter
plot(host..d15N~Light, data=D.win.df, col='blue', ylim=c(2, 8), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab="", cex.axis=0.8, cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1]),
           coef(summary(mod))["Light",1], 
      col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)

dev.copy(pdf, "figures/isotope/PanKB.d15N-host.pdf", width=6, height=4, encod="MacRoman")
dev.off()

Symbiont tissue nitrogen isotope composition

##################
# Fig: d15N symb over season and depth
##################        
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# model.data$symb..d15N
mod<-lmer(symb..d15N~Season+Light+Dom+ Season:Dom+ Light:Dom +(1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy_lmer(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params

# C summer
plot(symb..d15N~Light, data=C.sum.df, col='red3', ylim=c(2, 8), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon", 
           x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
par(new=T)

# D summer
plot(symb..d15N~Light, data=D.sum.df, col='blue', ylim=c(2, 8), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab=(expression(paste(delta^{15}, N[S], " (\u2030, Air)"))), cex.axis=0.8,
     cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]), 
           (coef(summary(mod))["Light",1]+ coef(summary(mod))["Light:DomD",1]), col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)

# C winter
plot(symb..d15N~Light, data=C.win.df, col='red3', ylim=c(2, 8), xlim=c(15, 0), yaxt="n", xaxt="n",
     cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]), coef(summary(mod))["Light",1], col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)

# D winter
plot(symb..d15N~Light, data=D.win.df, col='blue', ylim=c(2, 8), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab="", cex.axis=0.8, cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1]+
              coef(summary(mod))["Seasonwinter:DomD",1]),
           (coef(summary(mod))["Light",1]+ coef(summary(mod))["Light:DomD",1]), 
      col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)

dev.copy(pdf, "figures/isotope/PanKB.d15N-symb.pdf", width=6, height=4, encod="MacRoman")
dev.off()

Host-symbiont tissue nitrogen isotope composition

##################
# Fig: d15N host.symb over season and depth
##################        
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# model.data$d15N..host.sym
mod<-lmer(d15N..host.sym~Season+Light+Dom+ Season:Dom +(1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy_lmer(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params

# C summer
plot(d15N..host.sym~Light, data=C.sum.df, col='red3', ylim=c(-1,2), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon", 
           x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
par(new=T)

# D summer
plot(d15N..host.sym~Light, data=D.sum.df, col='blue', ylim=c(-1,2), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab=(expression(paste(delta^{15}, N[H-S], " (\u2030, Air)"))), cex.axis=0.8,
     cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]), 
           coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)

# C winter
plot(d15N..host.sym~Light, data=C.win.df, col='red3', ylim=c(-1,2), xlim=c(15, 0), yaxt="n", xaxt="n",
     cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]), coef(summary(mod))["Light",1], col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)

# D winter
plot(d15N..host.sym~Light, data=D.win.df, col='blue', ylim=c(-1,2), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), 
     ylab="", cex.axis=0.8, cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1]+
              coef(summary(mod))["Seasonwinter:DomD",1]),
           coef(summary(mod))["Light",1], 
      col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)

dev.copy(pdf, "figures/isotope/PanKB.d15Nh-s.pdf", width=6, height=4, encod="MacRoman")
dev.off()

Host tissue C:N

##################
# Fig: C.N host season and depth
##################       
      
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# log(model.data$host..C.N)
mod<-lmer(host..C.N~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy_lmer(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params

##### back transformed
x_light=seq(1,23, 0.5)
ypred= exp(intercept+ coef(summary(mod))["Light",1]*x_light)

# C summer
plot(host..C.N~Light, data=C.sum.df, col='red3', ylim=c(5, 12), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8,
     pch=21, bg="salmon", ylab="", xlab ="")
lines(x_light, ypred, col="salmon", lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
par(new=T)

# D summer
x_light=seq(5,24, 0.5)
ypred= exp(intercept+ coef(summary(mod))["DomD",1]+coef(summary(mod))["Light",1]*x_light)

plot(host..C.N~Light, data=D.sum.df, col='blue', ylim=c(5, 12), xlim=c(25, 0), cex=0.8, pch=21,
     bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
     ylab=(expression(paste("C:N"[H]), sep="")), 
     cex.axis=0.8, cex.lab=0.8, main="Summer", cex.main=0.7)
lines(x_light, ypred, col="deepskyblue1", lwd=2)


# C winter
x_light=seq(0.3,14, 0.5)
ypred= exp(intercept + coef(summary(mod))["Seasonwinter",1] + coef(summary(mod))["Light",1]*x_light)

plot(host..C.N~Light, data=C.win.df, col='red3', ylim=c(5, 12), xlim=c(15, 0), yaxt="n", xaxt="n",cex=0.8,
     pch=21, bg="salmon", ylab="", xlab ="")
lines(x_light, ypred, col="salmon", lwd=2)

par(new=T)

# D winter
x_light=seq(1.5,14, 0.5)
ypred= exp(intercept+ coef(summary(mod))["Seasonwinter",1] + coef(summary(mod))["DomD",1] +
             coef(summary(mod))["Light",1]*x_light)

plot(host..C.N~Light, data=D.win.df, col='blue', ylim=c(5, 12), xlim=c(15, 0), cex=0.8, pch=21,
     bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
     ylab="",
     cex.axis=0.8, cex.lab=0.8, main="Winter", cex.main=0.7)
lines(x_light, ypred, col="deepskyblue1", lwd=2)

dev.copy(pdf, "figures/isotope/PanKB.C.Nhost.pdf", width=6, height=4, encod="MacRoman")
dev.off()

Symbiont tissue C:N

##################
# Fig: C.N symb season and depth
##################       
      
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# log(model.data$host..C.N)
mod<-lmer(symb..C.N~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy_lmer(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params

##### back transformed
x_light=seq(1,23, 0.5)
ypred= exp(intercept+ coef(summary(mod))["Light",1]*x_light)

# C summer
plot(symb..C.N~Light, data=C.sum.df, col='red3', ylim=c(5, 12), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8,
     pch=21, bg="salmon", ylab="", xlab ="")
lines(x_light, ypred, col="salmon", lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
par(new=T)

# D summer
x_light=seq(5,24, 0.5)
ypred= exp(intercept+ coef(summary(mod))["DomD",1]+coef(summary(mod))["Light",1]*x_light)

plot(symb..C.N~Light, data=D.sum.df, col='blue', ylim=c(5, 12), xlim=c(25, 0), cex=0.8, pch=21,
     bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
     ylab=(expression(paste("C:N"[S]), sep="")), 
     cex.axis=0.8, cex.lab=0.8, main="Summer", cex.main=0.7)
lines(x_light, ypred, col="deepskyblue1", lwd=2)


# C winter
x_light=seq(0.3,14, 0.5)
ypred= exp(intercept + coef(summary(mod))["Seasonwinter",1] + coef(summary(mod))["Light",1]*x_light)

plot(symb..C.N~Light, data=C.win.df, col='red3', ylim=c(5, 12), xlim=c(15, 0), yaxt="n", xaxt="n",cex=0.8,
     pch=21, bg="salmon", ylab="", xlab ="")
lines(x_light, ypred, col="salmon", lwd=2)

par(new=T)

# D winter
x_light=seq(1.5,14, 0.5)
ypred= exp(intercept+ coef(summary(mod))["Seasonwinter",1] + coef(summary(mod))["DomD",1] +
             coef(summary(mod))["Light",1]*x_light)

plot(symb..C.N~Light, data=D.win.df, col='blue', ylim=c(5, 12), xlim=c(15, 0), cex=0.8, pch=21,
     bg="deepskyblue1",
     xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
     ylab="",
     cex.axis=0.8, cex.lab=0.8, main="Winter", cex.main=0.7)
lines(x_light, ypred, col="deepskyblue1", lwd=2)

dev.copy(pdf, "figures/isotope/PanKB.C.Nsym.pdf", width=6, height=4, encod="MacRoman")
dev.off()

Changes in host C and N isotope values show significant relationships, with significant relationships for C (summer and winter, p < 0.016) and D (summer, p =0.035) symbionts. No significant relationship for clade D in the winter was observed.

Changes in symbiont C and N isotope values were significant: in the winter for clade C (p=0.026) summer for clade D (p=0.021). Clade C in the summer and D in the winter were not significant.

Together this shows a significant relationship between carbon and nitrogen isotope values in clade D host and symbiont tissues during summer months, whereas the this relationship in clade C changes in both summer and winter (host only) and symbionts in the winter corals.

par(mfrow=c(1,4),mar=c(5,4.5,2,0.5))

summary(lm(host..d15N~host..d13C, data=C.sum.df)) # SIGNIF p=0.016
summary(lm(host..d15N~host..d13C, data=D.sum.df)) # SIGNIF p=0.035
summary(lm(host..d15N~host..d13C, data=C.win.df)) # SIGNIF p=0.004
summary(lm(host..d15N~host..d13C, data=D.win.df)) # not significant

# C summer plot
plot(host..d15N~host..d13C, data=C.sum.df, col='red3', xlim=c(-20, -12), ylim=c(2.5, 6.5), cex=0.7, pch=21, bg="salmon", main="Summer--host", cex.main=0.7, xaxt="n",
     ylab=(expression(paste(delta^{15}, N, " (\u2030, Air)"))),
     xlab=(expression(paste(delta^{13}, C, " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.9,
       y.intersp = 1.2, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.9, y.intersp = 1.2, x.intersp = 0.4, bty="n", inset=c(0.49, 0))
ablineclip(lm(host..d15N~host..d13C, data=C.sum.df), col='salmon', x1 = min(C.sum.df$host..d13C), x2 = max(C.sum.df$host..d13C), lwd=2)

# D summer plot
par(new=T)
plot(host..d15N~host..d13C, data=D.sum.df, col='blue', xlim=c(-20, -12), ylim=c(2.5, 6.5), cex=0.7, pch=21, bg="deepskyblue1", xaxt="n", yaxt="n", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(host..d15N~host..d13C, data=D.sum.df), col='deepskyblue1', x1 = min(D.sum.df$host..d13C), x2 = max(D.sum.df$host..d13C), lwd=2)

# C winter plot
plot(host..d15N~host..d13C, data=C.win.df, col='red3', xlim=c(-20, -12), ylim=c(2.5, 6.5), cex=0.7, pch=21, bg="salmon", xaxt="n", yaxt="n", main="Winter--host", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(host..d15N~host..d13C, data=C.win.df), col='salmon', x1 = min(C.win.df$host..d13C), x2 = max(C.win.df$host..d13C), lwd=2)

# D winter plot
par(new=T)
plot(host..d15N~host..d13C, data=D.win.df, col='blue', xlim=c(-20, -12), ylim=c(2.5, 6.5), cex=0.7, pch=21, bg="deepskyblue1", ylab="", cex.main=0.7, yaxt="n", xaxt="n",
xlab=(expression(paste(delta^{13}, C, " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
Axis(side=2, labels=FALSE)
ablineclip(lm(host..d15N~host..d13C, data=D.win.df), col='deepskyblue1', x1 = min(D.win.df$host..d13C), x2 = max(D.win.df$host..d13C), lwd=2)

#### r C and N isotope scatter symb

summary(lm(symb..d15N~symb..d13C, data=C.sum.df)) # not significant
summary(lm(symb..d15N~symb..d13C, data=D.sum.df)) # SIGNIF p=0.021
summary(lm(symb..d15N~symb..d13C, data=C.win.df)) # SIGNIF p=0.026
summary(lm(symb..d15N~symb..d13C, data=D.win.df)) # not significant

# C summer plot
plot(symb..d15N~symb..d13C, data=C.sum.df, col='palevioletred4', xlim=c(-20, -12), ylim=c(2.5, 6.5), cex=0.7, pch=21, bg="palevioletred2", main="Summer--symbiont", cex.main=0.7, ylab="", yaxt="n", xaxt="n",
     xlab=(expression(paste(delta^{13}, C, " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("palevioletred3", "lightskyblue3"),
       cex=0.9, y.intersp = 1.2, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("palevioletred2", "lightskyblue2"), col=c("palevioletred4", "lightskyblue4"), cex=0.9, y.intersp = 1.2, x.intersp = 0.4, bty="n", 
       inset=c(0.49, 0))
ablineclip(lm(symb..d15N~symb..d13C, data=C.sum.df), col='palevioletred2', x1 = min(C.sum.df$symb..d13C), x2 = max(C.sum.df$symb..d13C), lwd=2)

# D summer plot
par(new=T)
plot(symb..d15N~symb..d13C, data=D.sum.df, col='lightskyblue4', xlim=c(-20, -12), ylim=c(2.5, 6.5), cex=0.7, pch=21, bg="lightskyblue2", xaxt="n", yaxt="n", ylab="", xlab ="", cex.main=0.7)
Axis(side=2, labels=FALSE)
ablineclip(lm(symb..d15N~symb..d13C, data=D.sum.df), col='lightskyblue3', x1 = min(D.sum.df$symb..d13C), x2 = max(D.sum.df$symb..d13C), lwd=2)

# C winter plot
plot(symb..d15N~symb..d13C, data=C.win.df, col='palevioletred4', xlim=c(-20, -12), ylim=c(2.5, 6.5), cex=0.7, pch=21, bg="palevioletred2", xaxt="n", yaxt="n", ylab="", xlab ="", main="Winter--symbiont", cex.main=0.7)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
ablineclip(lm(symb..d15N~symb..d13C, data=C.win.df), col='palevioletred2', x1 = min(C.win.df$symb..d13C), x2 = max(C.win.df$symb..d13C), lwd=2)

# D winter plot
par(new=T)
plot(symb..d15N~symb..d13C, data=D.win.df, col='lightskyblue4', xlim=c(-20, -12), ylim=c(2.5, 6.5), cex=0.7, pch=21, bg="lightskyblue2", ylab="", xaxt="n", yaxt="n", cex.main=0.7,
xlab=(expression(paste(delta^{13}, C, " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=2, labels=FALSE)
ablineclip(lm(symb..d15N~symb..d13C, data=D.win.df), col='lightskyblue3', x1 = min(D.win.df$symb..d13C), x2 = max(D.win.df$symb..d13C), lwd=2)

dev.copy(pdf, "figures/regressions/PanKB.CN.scat.pdf", width=7, height=3, encod="MacRoman")
dev.off()
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

shallow.sum<-df.fig[df.fig$depth.bin.small=="<4m" & df.fig$Season=="summer", ]
deep.sum<-df.fig[df.fig$depth.bin.small==">4m" & df.fig$Season=="summer",]

shallow.wint<-df.fig[df.fig$depth.bin.small=="<4m" & df.fig$Season=="winter", ]
deep.wint<-df.fig[df.fig$depth.bin.small==">4m" & df.fig$Season=="winter",]


####### hosts
### shallow summer
plot(host..d15N~host..d13C, data=shallow.sum, col='red3', xlim=c(-20, -12), ylim=c(2.5, 6.5), cex=0.7, pch=21, bg="salmon", main="Summer", cex.main=0.7, xaxt="n",
     ylab=(expression(paste(delta^{15}, N[H], " (\u2030, Air)"))),
     xlab=(expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
legend("topright", c(" < 4m", " > 4m"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.31, 0))
ablineclip(lm(host..d15N~host..d13C, data=shallow.sum), col='salmon', x1 = min(shallow.sum$host..d13C), x2 = max(shallow.sum$host..d13C), lwd=2)

# deep summer
par(new=T)
plot(host..d15N~host..d13C, data=deep.sum, col='blue', xlim=c(-20, -12), ylim=c(2.5, 6.5), cex=0.7, pch=21, bg="deepskyblue1", xaxt="n", yaxt="n", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(host..d15N~host..d13C, data=deep.sum), col='deepskyblue1', x1 = min(deep.sum$host..d13C), x2 = max(deep.sum$host..d13C), lwd=2)


### shallow winter
plot(host..d15N~host..d13C, data=shallow.wint, col='red3', xlim=c(-20, -12), ylim=c(2.5, 6.5), cex=0.7, pch=21, bg="salmon", main="Winter", cex.main=0.7, xaxt="n", ylab="",
     xlab=(expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
legend("topright", c(" < 4m", " > 4m"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.31, 0))
ablineclip(lm(host..d15N~host..d13C, data=shallow.wint), col='salmon', x1 = min(shallow.wint$host..d13C), x2 = max(shallow.wint$host..d13C), lwd=2)

# deep winter
par(new=T)
plot(host..d15N~host..d13C, data=deep.wint, col='blue', xlim=c(-20, -12), ylim=c(2.5, 6.5), cex=0.7, pch=21, bg="deepskyblue1", xaxt="n", yaxt="n", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(host..d15N~host..d13C, data=deep.wint), col='deepskyblue1', x1 = min(deep.wint$host..d13C), x2 = max(deep.wint$host..d13C), lwd=2)

dev.copy(pdf, "figures/regressions/PanKB.scatdepth.host.pdf", width=6, height=4, encod="MacRoman")
dev.off()


#### symbionts

### shallow summer
plot(symb..d15N~symb..d13C, data=shallow.sum, col='red3', xlim=c(-20, -12), ylim=c(2.5, 6.5), cex=0.7, pch=21, bg="salmon", main="Summer", cex.main=0.7, xaxt="n",
     ylab=(expression(paste(delta^{15}, N[S], " (\u2030, Air)"))),
     xlab=(expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
legend("topright", c(" < 4m", " > 4m"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.31, 0))
ablineclip(lm(symb..d15N~symb..d13C, data=shallow.sum), col='salmon', x1 = min(shallow.sum$symb..d13C), x2 = max(shallow.sum$symb..d13C), lwd=2)

# deep summer
par(new=T)
plot(symb..d15N~symb..d13C, data=deep.sum, col='blue', xlim=c(-20, -12), ylim=c(2.5, 6.5), cex=0.7, pch=21, bg="deepskyblue1", xaxt="n", yaxt="n", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(symb..d15N~symb..d13C, data=deep.sum), col='deepskyblue1', x1 = min(deep.sum$symb..d13C), x2 = max(deep.sum$symb..d13C), lwd=2)


### shallow winter
plot(symb..d15N~symb..d13C, data=shallow.wint, col='red3', xlim=c(-20, -12), ylim=c(2.5, 6.5), cex=0.7, pch=21, bg="salmon", main="Winter", cex.main=0.7, xaxt="n", ylab="",
     xlab=(expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
legend("topright", c(" < 4m", "> 4m"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.31, 0))
ablineclip(lm(symb..d15N~symb..d13C, data=shallow.wint), col='salmon', x1 = min(shallow.wint$symb..d13C), x2 = max(shallow.wint$symb..d13C), lwd=2)

# deep winter
par(new=T)
plot(symb..d15N~symb..d13C, data=deep.wint, col='blue', xlim=c(-20, -12), ylim=c(2.5, 6.5), cex=0.7, pch=21, bg="deepskyblue1", xaxt="n", yaxt="n", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(symb..d15N~symb..d13C, data=deep.wint), col='deepskyblue1', x1 = min(deep.wint$symb..d13C), x2 = max(deep.wint$symb..d13C), lwd=2)

dev.copy(pdf, "figures/regressions/PanKB.scatdepth.symb.pdf", width=6, height=4, encod="MacRoman")
dev.off()

What is driving differences in d13C? Plot relationship between chlorophyll, chlorophyll per cell, and zoox in relation to d13C

  • In the host:zoox density in winter (C and D corals) positively relates with d13C values. In clade C corals, chl/cell influences d13C values in summer and shows near significant values in winter.

  • In the symbiont: we see similar patterns. In winter, d13C in clade C and D corals are influenced by zoox density. In both seasons for clade C corals, chl/cell also influences d13C values.

  • In host-symbiont: we also see that symbiont density relates to host-symb d13C. In winter, C and D corals show a relationship, but only D corals in the summer

### IN THE HOST ###

##### clade C
# ZOOX in the winter influence d13C
# CHLCELL in summer influences d13C

summary(lm(zoox~host..d13C, data=C.sum.df)) # NS
summary(lm(zoox~host..d13C, data=C.win.df)) # signif p=0.0001

summary(lm(chltot~host..d13C, data=C.sum.df)) # NS
summary(lm(chltot~host..d13C, data=C.win.df)) # NS

summary(lm(chlcell~host..d13C, data=C.sum.df)) # signif p = 0.024
summary(lm(chlcell~host..d13C, data=C.win.df)) # almost p=0.059


#### clade D
# ZOOX in winter influence d13C

summary(lm(zoox~host..d13C, data=D.sum.df)) # NS
summary(lm(zoox~host..d13C, data=D.win.df)) # signif p=0.030

summary(lm(chltot~host..d13C, data=D.sum.df)) # NS
summary(lm(chltot~host..d13C, data=D.win.df)) # NS

summary(lm(chlcell~host..d13C, data=D.sum.df)) # NS 
summary(lm(chlcell~host..d13C, data=D.win.df)) # NS


### IN THE SYMBIONT ###

##### clade C
# ZOOX in the winter influence d13C
# CHLCELL in summer and winter influences d13C

summary(lm(zoox~symb..d13C, data=C.sum.df)) # NS
summary(lm(zoox~symb..d13C, data=C.win.df)) # signif p<0.0001

summary(lm(chltot~symb..d13C, data=C.sum.df)) # NS
summary(lm(chltot~symb..d13C, data=C.win.df)) # NS

summary(lm(chlcell~symb..d13C, data=C.sum.df)) # signif p=0.004
summary(lm(chlcell~symb..d13C, data=C.win.df)) # signif p=0.030


#### clade D
# ZOOX in winter influence d13C

summary(lm(zoox~symb..d13C, data=D.win.df)) # signif p=0.007
summary(lm(zoox~symb..d13C, data=D.sum.df)) # NS

summary(lm(chltot~symb..d13C, data=D.win.df)) # NS
summary(lm(chltot~symb..d13C, data=D.sum.df)) # NS

summary(lm(chlcell~symb..d13C, data=D.win.df)) # NS
summary(lm(chlcell~symb..d13C, data=D.sum.df)) # NS 


##### IN HOST-SYMBIONT
## zoox
summary(lm(zoox~d13C..host.sym, data=C.sum.df)) # NS
summary(lm(zoox~d13C..host.sym, data=C.win.df)) # signif p=0.014
summary(lm(zoox~d13C..host.sym, data=D.sum.df)) # signif p=0.004
summary(lm(zoox~d13C..host.sym, data=D.win.df)) # signif p=0.007


############################
######### symbionts cm-2 and d13C

##### Plot it
##### host
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# C summer plot
plot(zoox~host..d13C, data=C.sum.df, col='red3', xlim=c(-20, -12), ylim=c(0, 6000000), cex=0.7, pch=21, bg="salmon", main="Summer", cex.main=0.7, xaxt="n",
     ylab=expression(paste("symbionts ", "(",cells~cm^-2~")", sep="")), 
     xlab=(expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
ablineclip(lm(zoox~host..d13C, data=C.sum.df), col='salmon', x1 = min(C.sum.df$host..d13C), x2 = max(C.sum.df$host..d13C), lwd=2)

# D summer plot
par(new=T)
plot(zoox~host..d13C, data=D.sum.df, col='blue', xlim=c(-20, -12), ylim=c(0, 6000000), cex=0.7, pch=21, bg="deepskyblue1", xaxt="n", yaxt="n", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(zoox~host..d13C, data=D.sum.df), col='deepskyblue1', x1 = min(D.sum.df$host..d13C), x2 = max(D.sum.df$host..d13C), lwd=2)

# C winter plot
plot(zoox~host..d13C, data=C.win.df, col='red3', xlim=c(-20, -12), ylim=c(0, 6000000), cex=0.7, pch=21, bg="salmon", xaxt="n", yaxt="n", main="Winter", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(zoox~host..d13C, data=C.win.df), col='salmon', x1 = min(C.win.df$host..d13C), x2 = max(C.win.df$host..d13C), lwd=2)

# D winter plot
par(new=T)
plot(zoox~host..d13C, data=D.win.df, col='blue', xlim=c(-20, -12), ylim=c(0, 6000000), cex=0.7, pch=21, bg="deepskyblue1", ylab="", cex.main=0.7, yaxt="n", xaxt="n",
xlab=(expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
Axis(side=2, labels=FALSE)
ablineclip(lm(zoox~host..d13C, data=D.win.df), col='deepskyblue1', x1 = min(D.win.df$host..d13C), x2 = max(D.win.df$host..d13C), lwd=2)

dev.copy(pdf, "figures/regressions/PanKB.zoox.d13Chost.pdf", width=6, height=4, encod="MacRoman")
dev.off()


####################
##### Plot it
##### symbiont
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# C summer plot
plot(zoox~symb..d13C, data=C.sum.df, col='red3', xlim=c(-20, -12), ylim=c(0, 6000000), cex=0.7, pch=21, bg="salmon", main="Summer", cex.main=0.7, xaxt="n",
     ylab=expression(paste("symbionts ", "(",cells~cm^-2~")", sep="")), 
     xlab=(expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
ablineclip(lm(zoox~symb..d13C, data=C.sum.df), col='salmon', x1 = min(C.sum.df$symb..d13C), x2 = max(C.sum.df$symb..d13C), lwd=2)

# D summer plot
par(new=T)
plot(zoox~symb..d13C, data=D.sum.df, col='blue', xlim=c(-20, -12), ylim=c(0, 6000000), cex=0.7, pch=21, bg="deepskyblue1", xaxt="n", yaxt="n", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(zoox~symb..d13C, data=D.sum.df), col='deepskyblue1', x1 = min(D.sum.df$symb..d13C), x2 = max(D.sum.df$symb..d13C), lwd=2)

# C winter plot
plot(zoox~symb..d13C, data=C.win.df, col='red3', xlim=c(-20, -12), ylim=c(0, 6000000), cex=0.7, pch=21, bg="salmon", xaxt="n", yaxt="n", main="Winter", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(zoox~symb..d13C, data=C.win.df), col='salmon', x1 = min(C.win.df$symb..d13C), x2 = max(C.win.df$symb..d13C), lwd=2)

# D winter plot
par(new=T)
plot(zoox~symb..d13C, data=D.win.df, col='blue', xlim=c(-20, -12), ylim=c(0, 6000000), cex=0.7, pch=21, bg="deepskyblue1", ylab="", cex.main=0.7, yaxt="n", xaxt="n",
xlab=(expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
Axis(side=2, labels=FALSE)
ablineclip(lm(zoox~symb..d13C, data=D.win.df), col='deepskyblue1', x1 = min(D.win.df$symb..d13C), x2 = max(D.win.df$symb..d13C), lwd=2)

dev.copy(pdf, "figures/regressions/PanKB.zoox.d13Csymb.pdf", width=6, height=4, encod="MacRoman")
dev.off()



############################
######### total chlorophyll a cm-2 and d13C

##### Plot it
##### host
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# C summer plot
plot(chltot~host..d13C, data=C.sum.df, col='red3', xlim=c(-20, -12), ylim=c(0, 12), cex=0.7, pch=21, bg="salmon", main="Summer", cex.main=0.7, xaxt="n",
     ylab=expression(paste("total chlorophyll", ~(mu*g~cm^-2), sep="")),
     xlab=(expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
ablineclip(lm(chltot~host..d13C, data=C.sum.df), col='salmon', x1 = min(C.sum.df$host..d13C), x2 = max(C.sum.df$host..d13C), lwd=2)

# D summer plot
par(new=T)
plot(chltot~host..d13C, data=D.sum.df, col='blue', xlim=c(-20, -12), ylim=c(0, 12), cex=0.7, pch=21, bg="deepskyblue1", xaxt="n", yaxt="n", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(chltot~host..d13C, data=D.sum.df), col='deepskyblue1', x1 = min(D.sum.df$host..d13C), x2 = max(D.sum.df$host..d13C), lwd=2)

# C winter plot
plot(chltot~host..d13C, data=C.win.df, col='red3', xlim=c(-20, -12), ylim=c(0, 12), cex=0.7, pch=21, bg="salmon", xaxt="n", yaxt="n", main="Winter", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(chltot~host..d13C, data=C.win.df), col='salmon', x1 = min(C.win.df$host..d13C), x2 = max(C.win.df$host..d13C), lwd=2)

# D winter plot
par(new=T)
plot(chltot~host..d13C, data=D.win.df, col='blue', xlim=c(-20, -12), ylim=c(0, 12), cex=0.7, pch=21, bg="deepskyblue1", ylab="", cex.main=0.7, yaxt="n", xaxt="n",
xlab=(expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
Axis(side=2, labels=FALSE)
ablineclip(lm(chltot~host..d13C, data=D.win.df), col='deepskyblue1', x1 = min(D.win.df$host..d13C), x2 = max(D.win.df$host..d13C), lwd=2)

dev.copy(pdf, "figures/regressions/PanKB.chltot.d13Chost.pdf", width=6, height=4, encod="MacRoman")
dev.off()


####################
##### Plot it
##### symbiont
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# C summer plot
plot(chltot~symb..d13C, data=C.sum.df, col='red3', xlim=c(-20, -12), ylim=c(0, 12), cex=0.7, pch=21, bg="salmon", main="Summer", cex.main=0.7, xaxt="n",
     ylab=expression(paste("total chlorophyll", ~(mu*g~cm^-2), sep="")),
     xlab=(expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
ablineclip(lm(chltot~symb..d13C, data=C.sum.df), col='salmon', x1 = min(C.sum.df$symb..d13C), x2 = max(C.sum.df$symb..d13C), lwd=2)

# D summer plot
par(new=T)
plot(chltot~symb..d13C, data=D.sum.df, col='blue', xlim=c(-20, -12), ylim=c(0, 12), cex=0.7, pch=21, bg="deepskyblue1", xaxt="n", yaxt="n", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(chltot~symb..d13C, data=D.sum.df), col='deepskyblue1', x1 = min(D.sum.df$symb..d13C), x2 = max(D.sum.df$symb..d13C), lwd=2)

# C winter plot
plot(chltot~symb..d13C, data=C.win.df, col='red3', xlim=c(-20, -12), ylim=c(0, 12), cex=0.7, pch=21, bg="salmon", xaxt="n", yaxt="n", main="Winter", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(chltot~symb..d13C, data=C.win.df), col='salmon', x1 = min(C.win.df$symb..d13C), x2 = max(C.win.df$symb..d13C), lwd=2)

# D winter plot
par(new=T)
plot(chltot~symb..d13C, data=D.win.df, col='blue', xlim=c(-20, -12), ylim=c(0, 12), cex=0.7, pch=21, bg="deepskyblue1", ylab="", cex.main=0.7, yaxt="n", xaxt="n",
xlab=(expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
Axis(side=2, labels=FALSE)
ablineclip(lm(chltot~symb..d13C, data=D.win.df), col='deepskyblue1', x1 = min(D.win.df$symb..d13C), x2 = max(D.win.df$symb..d13C), lwd=2)

dev.copy(pdf, "figures/regressions/PanKB.chltot.d13Csymb.pdf", width=6, height=4, encod="MacRoman")
dev.off()


############################
######### chlorophyll a cell-1 and d13C

##### Plot it
##### host
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# C summer plot
plot(chlcell~host..d13C, data=C.sum.df, col='red3', xlim=c(-20, -12), ylim=c(0, 12), cex=0.7, pch=21, bg="salmon", main="Summer", cex.main=0.7, xaxt="n",
     ylab=(expression(paste("chlorophyll (",pg~cell^-1, ")", sep=""))),
     xlab=(expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
ablineclip(lm(chlcell~host..d13C, data=C.sum.df), col='salmon', x1 = min(C.sum.df$host..d13C), x2 = max(C.sum.df$host..d13C), lwd=2)

# D summer plot
par(new=T)
plot(chlcell~host..d13C, data=D.sum.df, col='blue', xlim=c(-20, -12), ylim=c(0, 12), cex=0.7, pch=21, bg="deepskyblue1", xaxt="n", yaxt="n", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(chlcell~host..d13C, data=D.sum.df), col='deepskyblue1', x1 = min(D.sum.df$host..d13C), x2 = max(D.sum.df$host..d13C), lwd=2)

# C winter plot
plot(chlcell~host..d13C, data=C.win.df, col='red3', xlim=c(-20, -12), ylim=c(0, 12), cex=0.7, pch=21, bg="salmon", xaxt="n", yaxt="n", main="Winter", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(chlcell~host..d13C, data=C.win.df), col='salmon', x1 = min(C.win.df$host..d13C), x2 = max(C.win.df$host..d13C), lwd=2)

# D winter plot
par(new=T)
plot(chlcell~host..d13C, data=D.win.df, col='blue', xlim=c(-20, -12), ylim=c(0, 12), cex=0.7, pch=21, bg="deepskyblue1", ylab="", cex.main=0.7, yaxt="n", xaxt="n",
xlab=(expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
Axis(side=2, labels=FALSE)
ablineclip(lm(chlcell~host..d13C, data=D.win.df), col='deepskyblue1', x1 = min(D.win.df$host..d13C), x2 = max(D.win.df$host..d13C), lwd=2)

dev.copy(pdf, "figures/regressions/PanKB.chlcell.d13Chost.pdf", width=6, height=4, encod="MacRoman")
dev.off()


####################
##### Plot it
##### symbiont
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# C summer plot
plot(chlcell~symb..d13C, data=C.sum.df, col='red3', xlim=c(-20, -12), ylim=c(0, 12), cex=0.7, pch=21, bg="salmon", main="Summer", cex.main=0.7, xaxt="n",
     ylab=(expression(paste("chlorophyll (",pg~cell^-1, ")", sep=""))),
     xlab=(expression(paste(delta^{13}, C, " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.355, 0))
ablineclip(lm(chlcell~symb..d13C, data=C.sum.df), col='salmon', x1 = min(C.sum.df$symb..d13C), x2 = max(C.sum.df$symb..d13C), lwd=2)

# D summer plot
par(new=T)
plot(chlcell~symb..d13C, data=D.sum.df, col='blue', xlim=c(-20, -12), ylim=c(0, 12), cex=0.7, pch=21, bg="deepskyblue1", xaxt="n", yaxt="n", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(chlcell~symb..d13C, data=D.sum.df), col='deepskyblue1', x1 = min(D.sum.df$symb..d13C), x2 = max(D.sum.df$symb..d13C), lwd=2)

# C winter plot
plot(chlcell~symb..d13C, data=C.win.df, col='red3', xlim=c(-20, -12), ylim=c(0, 12), cex=0.7, pch=21, bg="salmon", xaxt="n", yaxt="n", main="Winter", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(chlcell~symb..d13C, data=C.win.df), col='salmon', x1 = min(C.win.df$symb..d13C), x2 = max(C.win.df$symb..d13C), lwd=2)

# D winter plot
par(new=T)
plot(chlcell~symb..d13C, data=D.win.df, col='blue', xlim=c(-20, -12), ylim=c(0, 12), cex=0.7, pch=21, bg="deepskyblue1", ylab="", cex.main=0.7, yaxt="n", xaxt="n",
xlab=(expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-20, -12, by = 2), cex.axis=0.8)
Axis(side=2, labels=FALSE)
ablineclip(lm(chlcell~symb..d13C, data=D.win.df), col='deepskyblue1', x1 = min(D.win.df$symb..d13C), x2 = max(D.win.df$symb..d13C), lwd=2)

dev.copy(pdf, "figures/regressions/PanKB.chlcell.d13Csymb.pdf", width=6, height=4, encod="MacRoman")
dev.off()

############################
######### zoox cm-2 and d13Ch-s

##### Plot it
##### host
par(mfrow=c(1,2),mar=c(5,4.5,2,2))

# C summer plot
plot(zoox~d13C..host.sym, data=C.sum.df, col='red3', xlim=c(-2, 2), ylim=c(0, 6000000), cex=0.7, pch=21, bg="salmon", main="Summer", cex.main=0.7, xaxt="n",
     ylab=expression(paste("symbionts ", "(",cells~cm^-2~")", sep="")), 
     xlab=(expression(paste(delta^{13}, C[H-S], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-2, 2, by = 1), cex.axis=0.8)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
       y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.35, 0))
ablineclip(lm(zoox~d13C..host.sym, data=C.sum.df), col='salmon', 
           x1 = min(C.sum.df$d13C..host.sym), x2 = max(C.sum.df$d13C..host.sym), lwd=2)

# D summer plot
par(new=T)
plot(zoox~d13C..host.sym, data=D.sum.df, col='blue', xlim=c(-2, 2), ylim=c(0, 6000000), cex=0.7, pch=21, bg="deepskyblue1", xaxt="n", yaxt="n", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(zoox~d13C..host.sym, data=D.sum.df), col='deepskyblue1', x1 = min(D.sum.df$d13C..host.sym), x2 = max(D.sum.df$d13C..host.sym), lwd=2)

# C winter plot
plot(zoox~d13C..host.sym, data=C.win.df, col='red3', xlim=c(-2, 2), ylim=c(0, 6000000), cex=0.7, pch=21, bg="salmon", xaxt="n", yaxt="n", main="Winter", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(zoox~d13C..host.sym, data=C.win.df), col='salmon', x1 = min(C.win.df$d13C..host.sym), x2 = max(C.win.df$d13C..host.sym), lwd=2)

# D winter plot
par(new=T)
plot(zoox~d13C..host.sym, data=D.win.df, col='blue', xlim=c(-2, 2), ylim=c(0, 6000000), cex=0.7, pch=21, bg="deepskyblue1", ylab="", cex.main=0.7, yaxt="n", xaxt="n",
xlab=(expression(paste(delta^{13}, C[H-S], " (\u2030, V-PDB)"))), cex.axis=0.8,
     cex.lab=0.8)
Axis(side=1, at = seq(-2, 2, by = 1), cex.axis=0.8)
Axis(side=2, labels=FALSE)
ablineclip(lm(zoox~d13C..host.sym, data=D.win.df), col='deepskyblue1', x1 = min(D.win.df$d13C..host.sym), x2 = max(D.win.df$d13C..host.sym), lwd=2)

dev.copy(pdf, "figures/regressions/PanKB.zoox.d13Chs.pdf", width=6, height=4, encod="MacRoman")
dev.off()
#Create a function to generate a continuous color palette
Pal <- colorRampPalette(c('gray80','gray30'))

#This adds a column of color values
# based on the y values

### one approach here
df.fig$Col <- Pal(10)[as.numeric(cut(df.fig$Light,breaks = 20))]

# plot(zoox~host..d13C, data=df.fig, pch = 20,col = df.fig$Col)


p1<-ggplot(df.fig, aes(x=host..d13C, y=zoox, color=newDepth)) + geom_point()+
scale_color_gradient(limits=c(0, 10), low="gray80", high="darkcyan") +
  ylab(expression(paste("symbionts", ~(cells~cm^-2), sep=""))) +
  xlab(expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)"))) +
  guides(color = guide_colourbar(barwidth = 0.8, barheight = 8)) + labs(color="Depth (m)")
  
p2<-ggplot(df.fig, aes(x=host..d13C, y=zoox, color=Light)) + geom_point()+
  scale_color_gradient(limits=c(0, 25), low="gray80", high="orange") +
  xlab(expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)"))) +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank())+
  guides(color = guide_colourbar(barwidth = 0.8, barheight = 8)) + labs(color="DLI")
  
grid.arrange(p1, p2, ncol = 2)

######################
#### zoox versus d13C-S

p3<-ggplot(df.fig, aes(x=symb..d13C, y=zoox, color=newDepth)) + geom_point()+
scale_color_gradient(limits=c(0, 10), low="gray80", high="darkcyan") +
  ylab(expression(paste("symbionts", ~(cells~cm^-2), sep=""))) +
  xlab(expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)"))) +
  guides(color = guide_colourbar(barwidth = 0.8, barheight = 8)) + labs(color="Depth (m)")
  
p4<-ggplot(df.fig, aes(x=symb..d13C, y=zoox, color=Light)) + geom_point()+
  scale_color_gradient(limits=c(0, 25), low="gray80", high="orange") +
  xlab(expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)"))) +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank()) +
  guides(color = guide_colourbar(barwidth = 0.8, barheight = 8)) + labs(color="DLI")
  
grid.arrange(p3, p4, ncol = 2)

######################
#### zoox versus d13C-HS

p5<-ggplot(df.fig, aes(x=d13C..host.sym, y=zoox, color=newDepth)) + geom_point()+
scale_color_gradient(limits=c(0, 10), low="gray80", high="darkcyan") +
  ylab(expression(paste("symbionts", ~(cells~cm^-2), sep=""))) +
  xlab(expression(paste(delta^{13}, C[H-S], " (\u2030, V-PDB)"))) +
  guides(color = guide_colourbar(barwidth = 0.8, barheight = 8)) + labs(color="Depth (m)")
  
p6<-ggplot(df.fig, aes(x=d13C..host.sym, y=zoox, color=Light)) + geom_point()+
  scale_color_gradient(limits=c(0, 25), low="gray80", high="orange") +
  xlab(expression(paste(delta^{13}, C[H-S], " (\u2030, V-PDB)"))) +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank())+
  guides(color = guide_colourbar(barwidth = 0.8, barheight = 8)) + labs(color="DLI")
  
grid.arrange(p5, p6, ncol = 2)

G<-arrangeGrob(p1,p2,p3,p4,p5,p6, ncol=2)
ggsave(file="figures/regressions/gradient.plots.png", G, height=8, width=8.5)